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on  the  idea  presented  in  (1),  of  having  a  two  Markov  chain  model 
to  describe  the  whole  system.  The  model  analyzed  in  this  paper  i 
simpler,  gives  further  insight  into  the  problem  of  interacting 
buffered  useis,  and  leads  to  one-dimentional  Markov  chains. 
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case  autocorrelation  curve  is  formulated. 
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An  infinite  cyclic  sequence  of  period  p,  where  p  S2 

for  some  positive  integer  k,  is  a  shift  register  sequence  if 

any  block  of  k  consecutive  bits  over  one  period  of  this 

k 

sequence  is  not  repeated.  The  sequence  of  period  p=2  -1  is 

known  as  the  maximal  length  sequence  and  the  sequence  of 

length  p=2k  is  called  a  deBruijn  sequence  [1].  For  a  given 

2k_1-k 

k  the  number  of  deBruijn  sequences  is  2 

A  linear  maximal  length  sequence  (also  known  as  a  PN 

sequence)  can  be  generated  by  a  k-stage  linear  feedback 

shift  register  (LFSR)  with  a  primitive  characteristic 

polynomial  [2],  For  a  given  k  the  number  of  linear  maximal 

u 

length  sequences  is  0(  2  -1)/k  .where  0  is  the  Euler's 
0-function.  An  efficient  method  to  construct  the  nonlinear 
sequences  is  the  cycle  joining  method  [2], [3]  which  is  based 
on  joining  the  sequences  of  short  period  together  to  obtain 
a  deBruijn  sequence. 

In  the  following  a  new  property  of  the  PN  sequences  is 
presented.  For  generation  of  nonlinear  sequences  of  period 
p=2  this  property  and  several  other  well  known  properties, 
[2],  of  the  PN  sequences  are  used  to  join  the  cycles  of  an 
LFSR  with  characteristic  polynomial  G ( x ) = g ( x ) ( 1 +x ) n  ,  where 
g(x)  is  a  primitive  polynomial  of  degree  m  and  G(x)  is  of 
degree  k=m+n. 


Super resolvi ng  image  restoration  (SIR)  in  the  presence  of 
noise  is  considered.  Few  '.*IR  algorithsm  have  demonstrated  the 
ability  to  resolve  two  point  sources  spaced  on-half  of  the 
Rayleigh  distance  apart.  In  this  report,  it  is  shown  that  the 
SIR  of  a  two-point  noncoherent  source  spaced  one-tenth  of  a 
Rayleigh  distance  apart  is  possible.  The  method  presented  uses 
optimal  data  fitting  techniques  based  on  the  method  of  linear 
programming.  For  noisy  images,  a  combination  of  linear 
eigenvalue  prefiltering  and  optimal  data  fitting  is  used.  It 
is  also  shown  that  for  a  d i f f ract ion-1 imi ted  (DL)  image  of  two 
point  sources  spaced  one-half  of  the  Rayleigh  distance  apart, 
where  the  input  is  contaminated  with  sigificant  noise,  SIR  is 
achievable.  These  results  have  important  implications  in 
atmospheric  physics,  geophysics,  radio  astronomy,  medical 
diagnostics,,  and  digital  bandwidth-compression  applications 
where  the  deconvolution  of  noisy  bandwidth-compressed  images 
is  one  of  the  fundamental  limitations. 

A  method  of  restoring  the  discrete  Fourier  transform  ( DFT) 
spectrum  of  a  DL  image  from  a  narrow  observation  segment  of  a 
DL  image  is  also  presented.  The  DL  spectral  restoration 
process  is  the  dual  of  the  more  common  DL  image  restoration 
process  with  the  role  of  frequency  and  space  reversed. 
Applications  of  a  spectrum  restoration  include  increasing  the 


field  of  view  of  existing  imaging  systems  and  extracting 
precise  frequency  components  of  a  large  DL  image  using  only 
small  segment  of  the  entire  image.  This  method  could  also  be 
employed  for  image  data  compression  which  is  of  interest  in 
digital  video  applications.  Several  differences  between  the 
implementation  of  the  image  and  the  spectrum  restoration 
processes  are  described.  The  estimate  is  constrained  to  have 
an  upper  bound  on  the  number  of  frequency  components  contained 
in  the  Fourier  spectrum.  The  bound  is  the  number  of  samples 
acquired  at  the  Nyquist  rate  for  the  length  of  the  image.  The 
magnitude  of  the  DFT  spectrum  is  also  bounded.  These 
constraints  define  a  large  number  of  possible  solutions.  The 
desired  solution  is  then  selected  such  that  the  distance, 
defined  in  a  function-theoretic  sense,  between  the  measured 
and  the  estimated  image  is  an  optimum.  A  number  such  measeures 
are  investigated.  Numerical  experiments  show  that  this 
approach  yields  results  that  are  highly  immune  to  measurement 
noise . 

Signal -dependent  noise  encountered  when  sampling  video  signals 
at  low  sampling  rate,  while  using  an  adaptive  video  delta 
modulator  both  as  a  source  encoder  and  a  video  digitizer,  is 
combated  by  using  a  white-light  reflective  transform  optical 
preprocessor.  While  the  white-light  preprocessor  works  on 
black  and  white  images,  its  main  advantage  is  that  it  can 
simultaneously  process  using  a  single  source  many  color 
channels.  Further,  the  relative  lack  of  both  spatial  and 


temporal  coherence  aids  in  reducing  speckle  and  interference 
noise  effects.  Experimental  results  on  color  video  images, 
that  have  been  preprocessed  using  a  white  light  optical 
transform  preprocessor  and  digitally  encoded  by  sampling  close 
to  the  Nyquist  rate  with  an  adaptive  delta  modulator,  are 
presented . 
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I.  Research  Objectives 


1.  Dynamic  Reservation  Multiple  Access  Technique 


for  Data  Transmission  via  Satellites 


A  new  reservation-based  mulitiple  access  packet  switch¬ 
ing  techniques  applicable  to  packet  communications  using  a 
satellite  channel  is  proposed.  The  objective  of  this  research 
is  to  optimize  the  system  performance  in  terms  of  throughput  and 
delay  and  to  develop  analytical  models  to  analyze  the  scheme. 


2 .  Analysis  of  Interacting  Buffered  Users  in  Slotted  ALOHA 


The  objective  of  this  research  is  to  analyze  the  slotted 
ALOHA  multiple  access  scheme  for  broadcase  channels  with  a  finite 
number  of  buffered  users.  Previous  investigators  of  this  problem 
have  assumed  that  a  user  has  a  maximum  of  one  packet  in  his  buffer. 
In  such  a  case,  if  a  second  packet  arrived  it  would  be  automatically 
rejected  and  destroyed. 


2) 


5.  Superresolving  Image  Restoration 

Superresolving  image  restoration  in  the  presence  of 
noise  is  considered. 

6 .  Restoration  of  the  DFT  Spectra 

Restoration  of  the  DFT  spectra  of  a  narrow  segment 
diffraction-limited  image  is  described. 

7 .  White-light  Prefiltering 
The  use  of  a  white-light  preprocessor  for  color  image 

i 

!  bandwidth  compression  is  described. 
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II.  Status  of  the  Research  Effort 
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DYNAMIC  RESERVATION  MULTIPLE  ACCESS  TECHNIQUE 
FOR  DATA  TRANSMISSION  VIA  SATELLITES 


D.  K.  Guha,  D.  L.  Schilling,  and  T.  N.  Saadawi 


Department  of  Electrical  Engineering,  City  College  of  New  York 
New  York,  New  York 


Abstract 


A  new  reservation- based  multiple  access  packet  switching  technique 
applicable  to  packet  communications  using  a  satellite  channel  is 
proposed  The  objective  of  this  research  is  to  optimize  the  system 
performance  in  terms  of  throughput  and  delay  aod  to  develop 
analytical  models  to  analyze  the  scheme 

Simulation  results  show  markedly  improved  dclay-ibroughput 
characteristics  over  other  existing  and  proposed  multiple  access 
techniques  with  especially  significant  performance  gains  at  high 
traffic  and  also  as  the  traffic  imbalance  among  the  users  increases. 
A  good  agreement  between  the  simulation  results  and  the  analytical 
results  is  obtained 

1.0  INTRODUCTION 

The  field  of  large  data  networks  has  seen  a  tremendous  growth 
in  the  last  decade.  The  simplest  solution  to  providing 
communication  between  two  points  is  to  assign  a  dedicated  channel 
for  their  use.  This  method  is  expensive  in  computer 
communications  especially  over  long  distances.  Measurement 
studies  conducted  on  lime-sharing  systems  indicate  that  both 
computer  and  terminal  data  streams  are  burny.  That  is,  the  peak 
data  rate  is  much  larger  than  the  average  data  rate.  Thus,  if  a 
channel  is  dynamically  shared  in  some  fashion  among  many  users, 
the  required  channel  capacity  may  be  much  less  than  the  unshared 
case  of  dedicated  channels.  This  concept  is  known  as  statistical  load 
averaging  and  has  been  applied  in  many  computer-  communication 
schemes  to  various  degrees  of  success  Tbe  application  of  packet 
switching  techniques  to  radio  communication  (both  satellite  and 
ground  radio  channels),  which  is  a  multi-access  broadcast  medium, 
provides  a  solution. 

In  the  subarea  of  multiple  access  schemes  for  networks  that 
include  broadcast  (and  particularly  satellite)  links  a  great  deal  of 
work  has  been  done  in  inventing  and  analyzing  strategics  and 
protocols  for  the  shared  use  of  a  common  channel  by  geographically 
separated  users.  The  problems  of  choosing  the  optimum  access 
protocol  (in  tbe  sense  of  minimizing  both  average  and  maximum 
delays  per  message  for  a  given  throughput)  have  not  been 
formulated  yet. 

The  basic  goal  of  this  paper  is  to  develop  control  schemes  that 
can  be  employed  when  using  a  geostationary  satellite  channel  for 
intercommunications  between  a  set  of  geographically  distributed 
nodes  and,  more  importantly,  modeling  and  analysis  of  such 
schemes  Although  many  lime-division  multiplexed  schemes  have 
been  proposed,  very  few  of  the  schemes  have  been  modeled 
analytically.  In  this  paper,  we  develop  and  analyze  a  dynamic 
reservation-based  multiple  access  packet  switching  technique,  which 
minimizes  both  average  and  maximum  delays  and  maximizes  traffic 
throughput.  The  analysis  is  hated  on  a  combination  of  two  Markov 


chain  model,  one  Markov  chain  describes  the  status  of  tbe  buffer 
contents  of  a  typical  user,  we  refer  to  it  as  the  User  Markov  Chain 
and  one  that  describes  the  tutus  of  all  the  users  of  the  channel,  we 
refer  to  it  as  the  System  Markov  Chain 

Simulation  resulu  have  been  obuined  for  the  system  under 
vanous  traffic  distributions  among  nodes,  and  are  presented  later  in 
the  paper  These  results  show  markedly  improved  delay- 
throughput  characteristics  over  other  existing  and  proposed  multiple 
access  techniques  with  especially  significant  performance  gains  as 
the  traffic  imbalance  among  the  users  increases.  A  good  agreement 
between  tbe  aimulauon  resulu  and  the  analytic  resulu  is  obuined 

For  equal  traffic  distribution  among  nodes,  each  node  needs 
approximately  two  (2)  slots  at  most  to  transmit  packeu  in  a  lime 
frame  to  achieve  the  best  performance.  It  is  found  that  the  average 
number  of  frame  size  is  approximately  I  slot  per  node  al  low 
throughput  and  less  than  2  slou  pet  node  at  high  throughput  In 
other  words,  the  knowledge  of  tbe  second  preceeding  reservation 
vector  for  a  node  is  not  necessary  to  have  two  slou  at  most  in  a 
time  frame 

For  unequal  traffic  distribution  among  nodes,  on  the  other  hand, 
the  delay-throughput  characteristic  it  improved  by  assigning 
variable  slots  to  the  nodes  in  a  time  frame.  The  knowledge  of 
preceding  reservation  vectors  is,  therefore,  necessary  to  assign 
variable  slou  to  a  node  in  a  time  frame  Optimum  values  of 
variable  slots  to  be  assigned  to  the  nodes  are  obuined  for  unequal 
traffic  among  nodes  in  which  tbe  arrival  rate  of  packeu  for  one 
node  is  eight  times  higher  than  that  of  other  nodes  However, 
these  values  will  depend  on  tbe  traffic  distribution  among  nodes 
Further  analysis  of  these  should  be  made. 

In  Section  2  0,  we  examine  several  existing  and  proposed 
satellite  communication  packet  switching  schemes.  A  new 
reservation  multiple  access  scheme  is  then  presented. 

In  Section  3.0.  this  proposed  scheme  is  modeled  and  analyzed  in 
a  queuing  theoretic  framework.  A  numerical  solution  to  the 
mathematical  model  is  obuined  Steady-Sute  System  performance 
parameters,  such  as  average  queue  sizes  and  waiting  times  which 
are  useful  in  system  design,  arc  obtained. 

In  Section  4.0,  a  FORTRAN  simulation  model  is  developed,  and 
a  comparison  is  made  between  the  analytical  solution  and  the 
simulated  performance  A  good  agreement  between  these  two  is 
obuined  The  simulation  resulu  are  compared  with  other  existing 
and  proposed  satellite  schemes  and  show  markedly  improved 
delay-throughput  characteristics  over  the  others.  Simulation  results 
arc  also  obuined  for  tbe  system  using  different  arrival  rales  of  two 
classes  of  traffic  under  vanous  loads  and  node  distributions.  These 
resulu  are  presented  later  in  this  section. 

In  Section  S.0,  conclusions  and  suggestions  for  further  work  are 
presented. 
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i.o  packet  switching  techniques  for  satelutfs 

1.1 

Several  techniques  have  been  previously  proposed  for  using  a 
satellite  channel  in  a  packet-switching  data  network  in  a  way  which 
allows  all  stations  to  dynamically  share  the  channel  capacity.  There 
is  a  variety  of  ways'  allocating  time  slots  to  the  nodes,  such  as 
fised-atsignmcnlt,  random  access,  reservation-techniques,  etc.  The 
well-known  ALOHA1  and  slolled-ALOHA1  achcmcs  are  eaamples 
of  random  access  techniques. 

Reservation-ALOHA4  used  random-accessing  in  conjunction 
with  a  frame  concept  and  dynamic  ownership  principle  to  allow  full 
utilization  of  the  channel  capacity.  A  different  reservation 
approach,  described  by  Roberts,’  uses  a  first-come  first-served 
distributed  queuing  system  in  conjunction  with  a  random-accessing 
technique  for  making  reservations. 

Binder's  technique*  consists  of  a  dynamic  time-assignment 
system  superimposed  on  a  fixed  STDM  structure.  The  approach 
retains  the  stability  and  allocation  fairness  of  a  filed  STDM  while 
allowing  nodes  to  make  use  of  channel  capacity  otherwise  wasted 
because  of  light  traffic  loads  and  a  non-uniform  distribution  of 
traffic  among  the  nodes. 

Balgangadhar  and  Pickholtz’  discuss  a  scheme  where  time  is 
divided  into  frames,  each  frame  consisting  of  reservation  time  slots, 
preassigned  time  slots  and  reicnrsiion  access  time  slots.  The  idea 
behind  this  scheme  is  to  make  efficient  channel  utilization  while 
circumventing  the  delay  problem  by  introducing  tbc  pre-assigned 
slots,  which  force  the  frame  size  to  be  no  smaller  than  one  round- 
trip  delay. 

1.1  A  Proposed  New  Reservation  Techalque 

The  proposed  scheme  consists  of  a  dynamic  reservation-based 
time-assignment  packet  twitching  technique.  In  this  proposed 
scheme,  lime  is  divided  into  frames,  each  frame  consisting  of 
reservation  time  slots,  preassigned  fixed  lime  slots  and  reservation- 
based  variable  time  slots.  At  the  beginning  of  each  frame,  there 
are  M  small  slots  (M  is  the  number  of  users  in  the  system),  each 
large  enough  to  transmit  a  reservation  from  one  of  the  M  nodes  in 
the  network  Following  these  small  reservation  slots,  there  are  M 
preassigned  fixed  slots,  with  each  node  being  permanently  allocated 
one  filed  slot  per  frame.  The  slot  size  and  the  number  of  nodes,  as 
in  Binder's  Scheme,  is  such  that  the  duration  of  these  M  fixed  slots 
is  at  least  as  long  at  one  round-trip  delay  to  the  satellite.  At  the 
end  of  these  fixed  slots  come  the  reservation  based  variable  slots 
(V).  The  number  of  these  variable  slots  in  each  frame  depends  on 
the  reservation  vectors  sent  by  the  nodes  at  the  beginning  of  tbe 
frames  and  tbe  algorithm  described  in  the  following  discussion. 
The  idea  behind  this  scheme  it  to  make  efficient  channel  utilization 
(hence,  the  asynchronous  variable  slots)  but  circumvent  the  delay 
problem  by  introducing  the  fixed  slots,  forcing  the  frame  size  to  be 
no  smaller  than  one  round-trip  delay.  One  major  advantage  of  this 
scheme  it  that  tbe  reservation  time  it  small  enough  to 
accommodate  two  bits  of  reservation  vector  for  each  node  and  thus 
a  very  small  percentage  of  the  total  frame  time,  which  it  not  true 
for  other  existing  schemes.  Figure  2.1  illustrates  Ihe  frame 
structure  and  time-slot  allocation. 

Tbe  reservation  vector  it  tent  by  each  node  at  follows:  The 
nodes  examine  Ibe  buffer  contents  (number  of  packets  opened)  just 
prior  to  sending  tbe  reservation  information,  subtract  one  from  it 
(since  they  know  that  they  will  get  one  fixed  slot  in  Ihe  frame)  and 
sends  two  bits  (C,0)  (or  only  0  in  this  case)  or  (0,1)  or  (1,0)  or 
(1,1)  into  the  two  reservation  slots  assigned  to  it.  Tbe  (0,0)  vector 
repreients  the  situation  where  the  station  does  not  have  any  packet 
for  transmission,  the  (0,1)  vector  represents  the  station  having  only 
one  packet  to  transmit,  the  (1,0)  vector  represents  the  station 
having  only  two  packets  to  transmit,  and  the  (I, I )  vector  represents 


the  station  having  more  than  two  packets  to  transmit.  Tbe  number 
of  variable  slots  assigned  to  each  station  in  a  time  frame  will  be 
determined  by  the  following  algorithm: 

A.  There  will  be  no  time  slot  allotted  for  the  station  which  has 
sent  a  (0,0)  vector  (or  only  ’0’  in  this  case). 

B  There  will  be  one  ume  slot  allotted  for  transmission  for  the 
station  which  has  sent  a  (0.1)  vector. 

C.  Thcrf  will  be  two  bme  slots  allotted  for  transmission  for  the 
station  which  has  sent  a  (1,0)  vector. 

D.  The  station  which  has  sent  a  (1,1)  vector  will  be  allotted 
time  slots  (equal  to  or  more  than  two)  which  will  depend  on 
its  prior  reservation  vectors  as  follows: 

I  P  (which  is  equal  to  or  more  than  two)  time  slots  will  be 
reserved  for  transmission  if  the  first  preceding 
reservation  vector  is  (0.0),  or  (0,1)  or  (1,0). 

2.  Q  (which  is  equal  to  or  more  than  P)  time  slots  will  be 
reserved  for  transmission  if  the  first  preceding  vector  is 
also  (1,1)  but  the  second  preceding  vector  is  (0.0)  or 
(0.1)  or  (1.0). 

3.  R  (which  is  equal  to  or  more  than  Q)  time  slots  will  be 
reserved  for  transmission  if  the  preceding  two  vectors 
arc  both  (1,1). 

Tbe  above  algorithm  is  explained  in  Table  2.1. 

Immediately,  after  all  the  nodes  have  transmitted  their 
reservation,  each  node  transmits  (if  it  has  any  packets  in  its  buffer), 
sequentially,  one  packet  on  Ibe  one  fised  slot  assigned  to  it  Sincaga 


Table  2.1 

Algorithm  for  Slot  Assignment 


RESERVATION  VECTORS  SENT  »T  EACH  N00E 

AT  THE  BEGINNING  Of  A  MIAMI 

NO  OF  RESERVATION 
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SLOTS  IV)  ASSIGNED 

TO  EACH  N00E  IN 

A  TIME  FRAME 
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Figure  2.1s.  Frame  structure  of  proposed  scheme 
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the  duration  of  the  M  fixed  slots  is  longer  than  a  round-trip  dciay, 
by  the  end  of  the  fixed  slots  all  the  nodes  will  have  received  the 
information  for  the  variable  slots  assigned  to  them.  Accordingly, 
each  node  will  transmit  its  assigned  packets,  sequentially,  in  the 
variable  portion  of  the  lime  frame. 

It  is  to  be  noted  that  the  slot  allocation  is  highly  asymmetrical 
with  regard  to  the  terminals.  It  it  intuitively  obvious  that  the 
packets  in  station  n  have  to  wait  longer  than  the  packets  of  station 
(n-l)  to  get  transmitted.  To  overcome  this  problem,  tbe 
assignment  of  the  slots  can  be  done  in  a  cyclic  order  in  frame  (j), 

tbe  sequence  of  slot  allocations  is  Q.2 . M);  in  frame  (j+1),  the 

sequence  is  changed  to  (2,3 . M.l).  and  to  forth;  in  frame 

(j+M+l),  the  sequence  it  again  (1,2 . M).  In  such  a  TDMA 

based  reservation  slot  tbe  reservation  minipacketa  are  quite  short, 
because  they  would  not  have  to  carry  identity  and  synchronization 
information. 

The  length  of  the  frame  is  a  random  variable  and  depends  on 
the  number  of  users  requesting  transmission  at  the  beginning  of  the 
frame.  The  mathematical  model  in  section  4  it  centered  primarily 
around  the  imbedded  Markov  chain,  imbedded  at  the  special 
epochs  which  arc  tbe  beginnings  of  the  frames.  The  service  order 
within  the  frame  has  no  effect  on  the  state  of  the  system  at  these 
epochs.  But  the  service  order  does  affect  the  waiting  time 
distribution  for  the  packets  queued  at  the  terminals. 

3.0  *  MATHEMATICAL  MODEL  OF  THE  PROPOSED  SCHEME 
3.1  Introduction 

The  proposed  system  is  a  non-sundard  multi-queuing  single¬ 
server  queuing  system.  Multiqueues  attended  by  a  single-server 
have  received  a  good  deal  of  attention  in  the  queuing  theory 
literature. 1-11  What  distinguishes  the  proposed  system  is  primarily 
the  service  discipline. 

Chu  and  Konbeim11  have  developed  a  unified  mathematical 
model  for  the  analysis  of  synchronous  TDM,  the  asynchronous 
TDM,  and  the  ’Hub  polling*  techniques.  Specifically,  they  assume 
that  the  arrival  pattern  of  packets  at  the  stations  it  Poisson  and, 
using  a  generating  function  approach,  have  developed  equations 
leading  to  the  queue  size  distributions  of  tbe  buffers  at  tbe  model 
and  the  waiting  times  experienced  by  the  packet.  , 

It  is  to  be  noted  that  tbe  service  time  of  a  packet  is  dependent 
on  the  status  of  the  remaining  queues.  If  we  let  nk  represent  the 
contents  of  the  input  buffer  of  the  kth  user,  then  a  complete 
description  of  all  M  users  requires  the  specification  of  tbe  joint 
probability  distribution  P(n|  ni,....nM)  The  determination  of  these 
probabilities  is  very  complex,  if  not  impossible  and  demands  the 
solution  of  a  large  number  of  seta  of  equations.  Fayolle  and 
lasnogorodski11  showed  the  limitation  c.  this  direct  approach.  They 
considered  the  simple  example  of  two  parallel  M/M/I  queuing 
systems  with  infinite  capacities  and  with  service  rate  for  each 
system  depending  on  the  status  of  the  other  system’s  queue  They 
showed  that  the  generating  function  F(x,  y),  corresponding  to  the 
joint  probability  of  the  two  queue  sizes  can  be  continued  as  a 
mcromorphic  function  to  the  whole  complex  plane.  Using  the 
theory  of  analytic  continuation  they  reduced  the  problem  to  a 
Rieman-Hilbert  problem  and  were  able  to  obtain  a  closed  form  for 
F(x,  y)  which  includes  several  elliptical  integrals  of  the  third  kind. 
Extension  of  their  analysis  to  the  case  of  more  than  two  users  and 
to  our  problem  seems  unfeasible. 

Leibowitz"  presents  an  approximate  method  of  treating 
multiqueuc  systems.  He  studied  tbe  case  of  N  queues  with  a  single 


server,  in  which  the  queues  arc  served  in  cyclic  order.  He  used  the 
following  argument  to  derived  the  probability  distribution  of  the 
queues;  if  the  server  meets  tbe  same  probability  distribution  of  the 
number  of  customers,  uy  P,.  ai  each  of  the  queues  on  one  cycle  of 
the  system,  then  the  same  distribution  must  meet  him  when  be 
returns  to  the  same  queue. 

Ftashida1’  used  Leibownz's  approximate  approach  in  the  analysis 
of  multiqueue  systems  where  the  server  serves,  at  mo  si.  K 
customers  that  were  waiting  when  he  arrived  at  a  queue  Also  in 
Reference  7,  the  authors  used  Leibowiu's  approximation  along 
with  an  independence  assumption  to  obtain  tbe  queue  sire 
distribution  for  each  user. 

Besides  introducing  a  new  reservation  scheme,  we  consider  that 
the  main  contribution  of  this  dissertation,  is  the  mathematical 
model  presented  here  to  analyze  interacting  buffered  terminals 
The  model,  as  mentioned  earlier,  is  basically  a  combination  of  two 
Markov  chains  imbedded  at  tbe  beginning  of  each  frame;  one 
Markov  chain  that  describes  the  cute  of  the  user,  referred  to  at  the 
User  Markov  Chain  and  another  Markov  chain  that  describes  the 
stale  of  all  the  users  in  the  systems,  referred  to  as  tbe  System 
Markov  Cham  T.  T.  Saadawi  and  A.  Ephremides1*  used  a  similar 
approach  to  analyze  the  reservation  scheme  where  they  considered 
that  each  station  transmits  maximum  one  packet  in  a  time  frame 

In  Section  3.2.1  we  describe  the  System  Markov  Chain  while  in 
Section  3.2.2  the  User  Markov  Chain  is  obtained.  For  the  purpose 
of  comparison  with  simulation  resulu,  the  mathematical  model  is 
studied  here  for  the  reservation  scheme  where  each  station  is 
allowed  to  transmit  maximum  two  packets  in  a  time  frame,  i.c.. 
P  —  Q  ••  R  - 1  For  other  values  of  P.  Q  and  R  we  will  have  similar 
analyses  for  the  system. 

3,2  Tbtl'sgr.Mgdtl 

The  system  consists  of  M  terminals,  each  of  which  has  an 
infinite  buffer.  The  arrival  process  at  each  of  the  M  terminals  is 
assumed  to  be  a  Bernoulli  process  with  rate  e,  i.e.,  tbe  probability 
of  arrival  of  a  (single  packet!  message  at  any  terminal  in  each  slot 
is  e  The  tout  arrival  rate  is,  therefore.  Me  packets  per  slot, 
which  is  equal  to  the  throughput  rale.  Most  other  studies  of 
multiple  access  protocols  have  assumed  a  Poisson  arrival  process 
The  Bernoulli  process  is  a  discrete-lime  analog  of  tbe  Poisson 
process,  which  is  well-suited  to  the  discrete  time  slot  structure. 
The  user  may  be  in  one  of  two  states,  an  idle  user,  where  bis  buffer 
is  empty,  or  an  active  user,  where  his  buffer  is  not  empty. 

As  shown  in  Fig  3  1  whenever  the  user  has  tbe  packet(s)  in  bis 
buffer,  it  sends  a  reservation  request  at  the  beginning  of  tbe 
frame*,  the  preassigned  packet  is  transmitted  to  the  fixed  time  slot 
and  the  reservation-based  packet(s)  are  transmitted  to  the  assigned 
slot(s)  of  the  variable  portion  of  the  frame  after  it  receives 
messages  from  the  satellite  at  tbe  end  of  a  round-trip  delay 

We  need  the  following  definitions; 

v,  —  steady  stale  probability  of  having  i  packets  in  the  buffer 
-  at  the  beginning  of  a  frame 

Hence 

-  probability  of  an  empty  buffer 

—  probability  of  an  idle  user 

I  —  T(  “  probability  of  an  active  user 


•  la  ih<  mathematical  model,  the  soistt  tesmeau  si  lbs  bepaaint  o{  *  frame, 
wbers  Hu  rttcrvsnoll  fooon  art  saal  hy  lbs  aodet,  sis  euamed  to  b*  of  sera 
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3.1.1  The  System  Markor  Chile  The  iUle  of  Ihc  system  it 
described  by  the  number  of  users  requesting  (nnsmiision  it  the 
beginning  of  the  frame.  Let  us  first  define  the  following: 

ji  —  number  of  users  hiving  only  one  picket  requesting 
transmission  it  the  beginning  of  i  frame 

j]  —  number  of  users  hiving  more  thin  one  picket 
requesting  transmission  it  the  beginning  of  i 
frame 

L>iri  “  length  of  i  frame  resulting  from  request!  for 
transmission  of  j,  ind  j,  users  it  its  beginning 

“  (M+jj)  slots,  where  M  is  cquil  to  or  grester  thin 
one  round-trip  deliy  in  time  slots 

"  steady  stite  probability  of  hiving  j,  and  j3  users 
requesting  transmission  it  the  beginning  of  a 
frame 


*i,~  probability  of  in  idle  user  generating  no  picket 
during  the  frame  of  length  L,^ 

-  (M0+i  (3.1.) 

“  probability  ef  in  idle  user  generating  one  picket 
during  the  frame  of  length  L^ 


M+j, 

I 


„  (l-n)^- 


(3.1b) 


h 


Fo- 


probability  of  in  idle  user  generating  more  thin 
me  picket  during  the  frame  of  length  L,^ 


(3.2) 


probability  of  hiving  two  pickets  only  in  the  buffer 
•t  the  beginning  of  i  frame  given  user  is  active  and 
hiving  more  thin  one  picket 


'i 


(3.3) 


F 


probability  of  hiving  three  pickets  in  the  buffer  it 
the  beginning  of  i  frame  giveo  user  is  active  and 
hiving  more  thin  two  pickets 


»i 


(34) 


In  the  above  binomial  distribution,  it  is  assumed  that  the 
probability  a  user  remains  active  is  independent  from  the  state  of 
the  other  users  A  similar  assumption  his  been  stated  and  justified 
in  Reference  15,  where  the  case  of  a  single  server  serving  N  lines 
in  cyclic  order  is  studied.  In  that  analysis,  the  number  of  lines  is 
issumed  constant  during  a  cycle  and  is  determined  by  a  binomial 
distribution. 


q(ii,i].Ki,K)J})  —  probability  |i,  idle  users  receive  one  packet 
and  i]  idle  users  receive  more  than  one 
packet  out  of  (M-K)  idle  users  during  a 
frame  given  K(—  K,-fK,)  users  remain  active 
out  of  j]  users  requesung  transmission  it  the 
beginning  of  i  frame] 
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Figure  3.1a.  System  Markov  Chain. 


The  System  Markov  Chain  is  illustrated  in  Fig.  3.1a.  Therefore, 
the  steady  stale  probability  of  n!  users  having  one  packet  and 
0]  users  having  more  than  one  picket,  requesting  transmission  it 
the  beginning  of  the  frame  is  given  by 

M  M'l| 

Fa, .a,  ”2  2  A„j,(|>i.,'j)P,la1  (3.7) 


where  Aj^fni.n,)  is  the  transition  probability  from  state  (ji jj)  to 
stile  (n,.n2)  and  is  giveo  by 


K|  —  number  of  users  having  one  packet  remain  active 
out  of  jj  users  requesting  transmission  it  the 
beginning  of  ■  frame 

Kj  —  number  of  users  hiving  more  than  one  picket 
remain  active  out  of  jj  users  requesting 
transmission  it  the  beginning  of  ■  frame 

M-K  —  number  of  idle  users  it  the  beginning  of  i  frame 
where  K  -  K,  +  K, 

C(K,.K]j])  —  probability  (K,  users  hiving  one  picket  ind  Kj 
users  hiving  more  thin  one  picket  remain  active 
out  of  j]  users  requesting  transmission  it  the 
beginning  of  a  frame) 


K,+K, 


(l-F.) 


fAN«, 


AWn»"*>  "2  2  MK  —  (Kj+Kj)  users  remain 

1,-0  Kp-0 

active  ,  ( n , — K , )  idle  users  receive  one  packet 
only,  and  (n2—  Kj)  idle  users  receive  more  than 
one  packet  during  a  frame,  given  j(—  ji+3i)  users 
request  transmission  at  the  beginning  of  a  frame.) 


2  2  q(n|— K„0j— Kj,K|,Kjj})C(Kl,KjjJ) 

>r°  a,-* 

for  1  £  j  —  j,  +  jj  S  M, 
and  Osn-ni  +  njSM 


for  j —  0,n2  «  0, 

ind  0£n~tt|  +  n}sM 
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3.1.1  The  Ultr  Markoi  Chain  The  stale  of  the  user  is  described 
by  the  number  of  packets  in  the  buffer  at  the  beginning  of  the 
frame. 

Let 

n  “  number  of  packets  in  the  buffer  at  the  beginning  of  a 
frame 

—  steady  sute  probability  of  having  n  packets  in  the  buffer 
at  the  beginning  of  a  frame 


•  a..ilnl 


('•I/ - •  •*•«»  t*> 


•  ) 


Figure  3.1b.  User  Markov  Chain. 

In  Fig.  3. 1.b,  we  show  the  User  Markov  Chain.  Note  that 
B„_i(n)  is  the  transition  probability  from  slate  (n-i)  to  state  (n). 
Let  us  now  determine  these  transition  probabilities.  Let 

g<  i.Lj|  jj)  —  Prti  packets  arrive  during  the  frame  of  length  Lj 


M+jj]  . 

-[  .  Ja*(l— a)"41’ 


B,4!( n)  —  Prlrwo  packets  leave  and  no  packet  arrives  during  the 
frame  of  length  L^^] 

M  M*il  I  M 

-2  2«<°.IW«W  l-2P,,.o.  n- 0.1.2... (3.10) 
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Similarly. 

Bt«i(n)  —  Prltwo  packets  leave  and  one  packet  arrives  during  the 
frame  of  length  Lj^l 

M  M-j,  {  M  \ 
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frame  of  length  L^j,) 

B,(0)  -  Prfone  packet  leaves  and  no  packet  arrives  during  the 
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Similarly, 


B.(n)  —  Prltwo  packet*  leave  and  two  packet*  arrive  during  tbe 
frame  of  length 

-  2  2«(2.^J,)P<(J/  1-2 p, 4  «  -  2.3.  ...  (3.13) 
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B,(l)  —  Prlonc  packet  leaves  and  one  packet  arrives  during  the 

frame  of  length 
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Bq(0)  -  Prlno  packet  arrives  during  the  frame  of  length  of  L,(Jj) 

-  “2  M2”l '-2Pm-,J  (3.35 
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Similarly, 

Br^(n)  -  Prltwo  packets  leave  and  (i+2)  packets  arrive  dunng 
the  frame  of  length  L^j,) 


-2  2 
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B|( n)  —  Prln  packets  arrive  during  the  frame  of  length  LMj] 
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Bg(n)  —  Pr  In  packets  arrive  during  the  frame  of  length 
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-  2  2  *(».Wpm/  '-2p«-,J-  (J  '*) 

h~°  h~*  {  *i"°  J 

n-  1 . 2M-1 


We  now  can  have  a  numerical  solution  for  »,.n  —  0.1.2,  3,...  as 
described  below. 

3.1.3  Numerical  Solstices  Equations  (3.1)  through  (3.6)  provide 
the  values  of  the  quantities  that  were  needed  in  order  to  solve  the 
equations  (3.7).  However,  some  of  these  quantities  arc  not 
eapressed  in  terms  of  only  the  system  parameter  Me.  but,  in  terms 
of  F8  and  F,  which  are  function  of  »i,  ejand  vj,  which  in  turn 
are  eapressed  in  terms  of  |PB|  a)l-  Thus,  in  total,  we  have  a  set  of 
simultaneous,  coupled  nonlinear  equations  in  r,  and  lP,t.,,l.  In 
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Fig  3.2.  «e  show  the  flow  chart  lo  solve  for  these  values.  We  iiart 
»nb  an  initial  value  for  F©  and  F,.  then  solve  equation  (3.7)  m 
Viing  the  values  of  Pl|(lj  we  then  determine  the  steady 
state  probabilities  r,  which  satisfy  the  following  maths  equation: 

-*.B  (3.19) 

where  B  is  the  transition  probability  matrix  with  elements 
b„-B,(j) 


The  average  frame  sue.  another  parameter  of  our  interest,  is 
given  by 


M  M“l| 


s<-2  2 

if*  if* 


(123) 


Finally.  we  define  aubiiity  as.  *T he  tysic  m  it  in  id  to  be  mbit  if 
the  iteady  state  prohnbiliiy  of  having  »  finuc  number  of  pnekeu  in 
Ibe  buffer  emu  ind  it  nonzero* 
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Figure  3.2.  Flow  chart  to  determine  average  delays  for 
the  proposed  scheme. 

The  solution  of  n  equations  in  n  unknowns  in  equation  (3.19)  is 
straightforward  using  matrix  inversion*  and  equations  (3.9)  through 
(3.18).  Then,  using  Wegestein's  iteration  scheme17,  we  obtain  the 
new  F0  and  Fi  and  repeat  it  until  it  converges  to  a  solution  for 
A  FORTRAN  program  was  written  to  solve  r.  The  analytical 
results  are  compared  with  the  simulation  result!  in  the  nest 
chapter. 

3.2.4  Delay  Analysis  The  average  total  delay  per  picket,  D,  is  the 
sum  of  the  average  waiting  time  in  the  buffer,  W,  plus  the  average 
time  from  the  beginning  of  the  frame  till  the  packet  ii  irammitted 
in  iu  assigned  slot.  We  refer  to  the  latter  as  the  service  time,  S. 
Hence, 

D  -  W  +  S  (3.20) 

Using  Little's  Result,"  we  can  write  W  u, 

W--2-  (3.21) 

where. 


«.»  SIMULATION  RESULTS 

A  FORTRAN  simulation  program  was  wntlcn  to  find  the 
throughput-delay  performance  of  the  proposed  model  described  in 
the  last  chapter,  under  the  following  traffic  conditions. 

a  Equal  traffic  distribution  among  nodes  in  which  traffic  was 
equally  distributed  among  the  nodes 
•  Unequal  traffic  distribution  among  nodes  in  which  the 
arrival  rate  of  packets  for  one  *big‘  node  is  higher  than  that 
of  other  'small'  nodes 

It  will  be  assumed  that  the  satellite  channel  is  a  highly  reliable  data 
link,  in  which  noise  may  be  neglected  (a  typical  value  for  the 
capacity  of  the  channel  it  $0,000  bits/scc).  Packets  of  fiaed  length 
are  used  and  lime  is  slotted  so  that  one  packet  is  transmitted  per 
slot.  The  slot  tire  and  the  number  of  nodes  is  such  that  the 
duration  of  these  M  fixed  slots  is  one  round-trip  delay  lo  the 
satellite  i.e..  .27  sec. 

4.1  Equal  Traffic  Amomt  NWfl 

In  this  case,  the  arrival  process  at  each  of  the  M  terminals  is 
assumed  to  be  a  Bernoulli  process  with  rate  a,  i.e..  the  probability 
of  arrival  of  a  (single  packet)  message  at  any  terminal  in  each  slot 
is  a.  The  total  arrival  rate  is.  therefore,  M  a  packets  per  slot, 
which  is  same  as  the  tysicm  throughput  or  Ibe  system  utilization 
factor.  The  Burnoulli  process  it  well  suited  to  the  discrete  time  slot 
structure  and  is  a  discrete-time  analog  of  the  Poisson  process.  For 
varying  M.  the  number  of  nodes  and  Mo,  the  system  utilization 
factor,  the  following  statistics  have  been  collected  from  the 
simulation  model.  In  a  finite  run-length  simulation,  the  statistics 
are  an  arithmetic  average  over  the  M  queues.  Also,  a  comparison 
between  these  system  performance  parameters  obtained  by  the 
mathematical  model  and  those  obtained  by  the  simulation  model 
hxi  been  made  for  M  —  4  for  the  xituation  where  each  node  is 
allowed  to  transmit  maximum  two  packets  in  a  time  frame  (i.e., 
P  -Q  -  R  -  1)  for 

(1)  Mean  and  standard  deviation  of  the  queue  sizes  at  the 
nodes  at  the  beginning  of  a  frame  (Table  4.1). 

(2)  Mean  and  standard  deviation  of  the  frame  length  expressed 
in  number  of  time  slots  (Table  4.2) 

(3)  Average  and  maximum  delay  for  a  packet  at  various  loads 
(Figures  4.2  and  4.3).  The  analytical  result  for  the  average 
delay  for  M“4  and  P—Q—R  —  I  is  also  compared  with  the 
simulation  result  in  Figure  4.2. 


Q  “  Average  Queue  Size 

«D 

-  2n  ‘  *• 
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Assuming  cyclic  assignment  discussed  earlier  in  Section  2.2,  the 
average  service  time  ii  given  by 
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*  When  •  u  very  lar|e.  il  alight  be  Silticuli  is  invert  Ike  aitrii  ■  la  that  caie, 
aiher  Uaa4ar4  methods  thaa!4  ha  lied  ta  salve  aqaaliaa  0.19) 


Figure  4.1  compares  the  average  delay  (or  the  proposed  scheme 
with  that  for  other  protocols  discussed  in  Section  2.0.  It  shows  that 
the  proposed  scheme  gives  a  significant  improvement  (or  delay  over 
the  other  reservation  schemes.  It  it  also  teen  from  Figures  (4.2 
and  4.3)  that,  for  the  traffic  situation  discussed  above,  each  node 
needs  approximately  two  (2)  tlott  at  most  to  transmit  packets  in  a 
time  frame  to  achieve  the  best  performance.  This  result  also  agrees 
with  the  results  of  Balgangadhar  and  Pickholtz's  tchcme7  where 
they  have  found  out  that  the  average  number  of  frame  size  is 
approximately  I  slot  per  node  at  low  throughput  and  less  than  2 
slots  per  node  at  high  throughput  (See  Table  4.2).  In  other  words, 


the  knowledge  of  the  lecond  preceding  reiervation  rector  for  ■ 
node  is  not  necessary  to  have  two  slots  at  most  in  a  time  frame. 

The  agreement  between  the  analytical  results  and  the  simulation 
results  is  seen  to  be  eaccllcnt  over  a  wide  nngc  of  Me,  the  system 
utilization  factor.  For  eaample,  for  M  •  4  and  the  utilization 
factor  of  0.4,  the  mean  queue  size  at  the  epochs  and  the  mean 
frame  size,  computed  analytically,  are  0.4]  and  4.24  rrraectively 
(see  Tables  4.1  and  4.2).  For  the  same  value  of  Me,  simulation 
results  give  the  mean  queue  size  at  the  epoch  as  0.44  and  the  mean 
queue  size  at  the  epochs  as  0.44  and  the  mean  frame  size  as  4.29. 
The  disagreement  between  these  and  the  analytical  values  it  about 
2  percent  and  1  percent,  respectively. 

The  agreement  between  the  analytical  and  the  simulation  results 
(Figure  4.2)  it  teen  to  be  particularly  good  for  low  values  of  tysttm 
utilization  factor.  It  is  difficult  to  argue  why  the  independence 
assumption  made  in  the  mathematical  analysis  should  hold  for  low 
value:  of  the  utilization  factor.  On  the  other  hand,  for  large  values 
of  M,  one  might  argue  that  because  of  the  large  number  of  queues, 
the  different  queues  become  decoupled  and  hence  make  the 
independence  assumption  valid.  But  further  analysis  of  this  should 
be  made. 

Table  4.1 

Comparison  Of  Computed  (Analytical)  and  Simulated  Values  for 
the  Mean  and  Standard  Deviation  (In  Parenthesis,  Below  the 
Mean)  of  Queue  Sizes  at  the  Beginning  of  a  Frame,  for  Various 
Values  of  N  and  Na  (the  Number  of  Nodes  and  the  System 
Utilization  Factor,  Respectively) 


Ns  ■  System  Utilization 


0.1 

0  2 

0.3 

04 

0  5 

0  6 

0  7 

0.8 

Aru1)iical 

I 

.2 

.31 

43 

.57 

.72 

.91 

1.15 

<M  -  4) 

(3) 

1.45) 

(54) 

(.43) 

(.7) 

(.77) 

(.14) 

(0.88) 

Stmubied 
M  -  4 

.1 

2 

32 

.44 

.57 

.72 

.93 

1.2 

(3) 

(.45) 

(.54) 

(■64) 

(72) 

(79) 

(  9) 

(M) 

M  -  12 

.1 

.2 

.31 

43 

54 

.72 

.91 

1.17 

(3) 

(.45) 

(.55) 

(.651 

(.75) 

(85) 

(96) 

(1.09) 

M  -  20 

.1 

.2 

.31 

42 

.55 

.74 

.92 

1.21 

(.31) 

(44) 

(.541 

(65) 

(.74) 

(85) 

(.97) 

(1.12) 

Figure  4.1.  Comparison  of  various  schemes  (12  nodes): 
average  packet  delay  vs.  throughput. 


Figure  4.2.  Average  Packet  Delay  vs.  Throughout  (4  Nodes) 


Table  4.2 

Comparison  Of  Computed  (Analytical)  and  Simulated  Values  for 
the  Mean  and  Standard  Deviation  (In  Parenthesis,  Below  the 
Mean)  of  Frame  Size  in  Number  of  Slots,  for  Various  Values  of  N 
and  N#  (the  Number  of  Nodes  and  the  System  Utilization  Factor, 
Respectively) 
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Figure  4.].  Maaimum  packet  delay  vs.  throughput  (4 
nodes). 
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4.1  Unequal  Traffic  Aphi  NWw 

In  this  case,  it  is  considered  lhal  a  ’big*  node  consists  of  (single 
packet)  messages  whose  arrival  rate  («,)  is  eight  limes  the  arrival 
rate  (a,)  of  (single  packet)  messages  at  other  ‘small*  nodes.  Total 
small  node  traffic  was  divided  equally  among  themselves.  The  total 
arrival  rate  it,  therefore,  (#k+(M-l)e,).  which  it  equal  to  the 
system  throughput. 

For  varying  M,  the  number  of  nodes,  and  the  system 
throughput,  the  average  (Figure  4.4)  and  maaimum  (Figure  4.S) 
delay  for  a  packet  at  'small*  and  ’big*  nodes  at  various  loads  have 
been  obtained  from  the  simulation  model.  It  it  teen  from  these 
figures  lhal  by  assigning  variable  slots  to  the  ‘big*  node  we  improve 
the  average  and  maaimum  delay  characteristics  of  the  system.  The 
knowledge  of  preceding  reservation  vectors,  it,  therefore,  necessary 
to  assign  variable  slots  to  a  node  in  a  time  frame. 


Figure  4.4.  Average  packet  delay  vs.  throughput  (4 
nodes)  packet  arrival  rate  of  I  big  node  eight 
times  higher  than  that  of  each  of  3  small 
nodes. 

The  same  algorithm  for  reservation  of  slots  is  applied  to  all 
nodes  (Table  3.1).  It  it  teen  that  the  best  performance  is  obtained 
for  P“3,  Q- 6,  and  R~ll  for  the  traffic  situation  discussed  above. 
However,  optimum  values  of  P,  Q  and  R  and  their  relationships 
will  depend  on  the  ratio  of  the  arrival  rates  at  ’big*  and  ‘small* 
nodes.  Further  analysis  of  these  should  be  made.  It  it  important 
to  note  here  that  an  adaptive  algorithm  (similar  to  the  algorithm 
discussed  in  this  dissertation),  where  a  node  asking  for  more  than 
two  slots  Iby  tending  reservation  vector  (1.1)1  is  assigned  some 
additional  slots  in  addition  to  the  slots  assigned  to  it  in  the  previous 
frame,  is  alto  studied  and  found  not  auitable  in  this  cue. 

5.0  CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  WORK 

The  field  of  large  data  networks  hat  teen  a  tremendous  growth 
in  the  lut  decade.  In  the  subarea  of  multiple  access  schemes  for 
neworks  that  include  broadcast  (and  particularly  satellite)  links, 
many  time-division  mulliplcaed  schemes  have  been  proposed  for 
the  shared  use  of  a  common  channel  by  geographically  separated 


Figure  4.5.  Maximum  packet  delay  vs.  throughput  (4 
nodes)  packet  arrival  rate  of  I  big  node  eight 
times  higher  than  that  of  each  of  3  small 
nodes 


users,  but  very  few  of  the  schemes  have  been  modeled  analytically. 
Also,  the  problems  of  choosing  tbc  optimum  access  protocol  (in  the 
tense  of  minimizing  average  delay  per  message  for  a  given 
throughput)  have  not  been  formulated  yet.  In  tbit  paper  we 
develop  optimum  control  ’  schemes  that  can  be  employed  when 
using  a  geostationary  satellite  channel  for  intercommunications 
between  a  set  of  geographically  distributed  nodes  and,  more 
importantly,  the  modeling  and  analysis  of  such  schemes 

Several  existing  and  proposed  satellite  multiple  access 
techniques  are  examined.  In  this  dissertation,  we  develop  and 
analyze  a  dynamic  reservation-based  multiple  access  packet 
switching  technique,  which  minimizes  both  average  and  maximum 
delays  and  maximizes  traffic  throughput.  The  analysis  is  taxed  on  a 
combination  of  two  Markov  chain  model;  one  Markov  chain 
describes  the  status  of  the  bufTer  contents  of  a  typical  user,  we  refer 
to  it  as  the  User  Markov  Chain  and  one  that  describes  the  status  of 
all  the  users  of  the  channel,  we  refer  to  it  as  tbc  System  Markov 
Chain. 

Simulation  results  have  been  obtained  under  various  traffic 
distributions  among  nodes.  These  results  show  markedly  improved 
delay-throughput  characteristics  over  other  existing  and  proposed 
multiple  access  techniques  with  especially  significant  performance 
gains  at  high  traffic  and  alto  as  tbc  traffic  imbalance  among  the 
users  increases.  A  good  agreement  between  the  simulation  results 
and  tbe  analytical  results  it  obtained. 

For  equal  traffic  distribution  among  nodes,  each  node  needs 
approximately  two  (2)  slots  at  most  to  transmit  packets  in  a  time 
frame  to  achieve  the  best  performance.  It  it  found  that  tbe  average 
number  of  frame  size  is  approximately  I  slot  per  node  at  low 
throughput  and  lets  than  2  slots  per  node  at  high  throughput.  In 
other  words,  the  knowledge  of  the  second  preceding  reservation 
vector  for  a  node  it  not  necessary  to  assign  two  slots  at  most  in  a 
time  frame. 


12) 


For  unequal  traffic  distribution  among  nodes,  on  tbe  other  hand, 
the  delay-throughput  characteristic  is  improved  by  assigning 
variable  slots  to  the  nodes,  in  a  lime  frame.  Tbe  knowledge  of 
preceding  reservation  vectors  is,  therefore,  accessary  to  assign 
variable  slots  to  a  node  in  a  lime  frame.  The  optimum  values  of 
variable  slots  to  be  assigned  to  the  nodes  are  obtained  for  unequal 
traffic  among  nodes  in  which  the  arrival  rate  of  packets  for  one 
node  is  eight  times  higher  than  that  of  other  nodes.  However, 


these  values  will  depend  on  the  traffic  distribution  among  nodes 
Further  analysis  of  these  should  be  made 

The  proposed  packet  switching  technique  based  upon  dynamic 
reservation  multiple  access  concept  may  also  be  applied  to  the  use 
of  ground  radio  channels  for  terminal-computer  communications 
In  Ihis  case  the  analysis  will  be  the  same,  except  the  time  frame 
will  consist  of  only  reservation-based  variable  time  slots  but  no 
preassigned  hied  slots  because  there  is  no  round-tnp  delay. 
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ABSTRACT  ‘L 

This  paper  analyzes  the  slotted-ALOHA  multiple  access  scheme 
for  broadcast  channels  with  a  finite  number  of  users,  each 
having  a  buffer  of  infinite  capacity.  The  analysis  is 
based  on  the  idea  presented  in  (1),  of  having  a  two  Markov 
chain  model  to  describe  the  whole  system.  The  model 
analyzed  in  this  paper  is  simpler,  gives  further  insight 
into  the  problem  of  interacting  buffered  users,  and  leads  to 
one-dimensional  Markov  chains. 

1 .  INTRODUCTION 

The  problem  of  analyzing  random  access  schemes  for  a 
finite  number  of  buffered  users  is  still  an  open  problem. 

In  most  of  the  previous  analysis  for  random  access 
schemes,  (2),  (3),  I1!)  and  (5),  it  was  assumed  that  a  user 

can  have  a  maximum  of  one  packet  in  his  buffer.  Any  packet 
that  is  generated,  while  there  is  a  packet  in  the  buffer,  if 
assumed  to  be  rejected. 

Buffering  is  essential  in  most  practical 
applications,  such  as  satellite  access  networks.  When 


buffering  is  not  used,  packet  rejection  is  often  employed. 
Packet  rejection  is  costly  and  represents  a  drastic  means  of 


flow  control  . 


In  this  paper,  and  as  a  typical  example  of  a  random 
access  scheme,  we  analyze  the  slotted  ALOHA  system  with  a 
finite  number  of  buffered  users.  The  approach  here  is 
similar  to  the  one  presented  in  (1),  yet  the  model  presented 
here  is  simpler  and  Rives  further  insip. ht  into  the  problem. 
The  approach  is  based  on  modelling  the  interaction  between 
the  competing  terminals  by  a  combination  of  two  Markov 
chains.  The  key  idea  behind  the  interaction  is  the  fact 
that  successful  transmission  of  a  packet  by  a  user  depends 
on  the  status  of  the  other  users.  Thus  it  is  possible  to 
consider  the  state  of  a  user  (represen ter)  by  the  size  •  :  ms 
queue)  modeled  by  an  imbedded  Markov  chain,  whose  p  a  r  a  ••  r. 

(transition  probabilities)  depend  on  t.lv  state  <  •  •  .••••- 

(basically,  the  number  of  user.-,  who  . . . 

those  who  do  not).  The  latter  state  . 

represented  by  an  imbedded  Markov  chain  wt..,se  ,  ■  f  i  - . --  •  • 
in  turn  functions  of  the  state  of  the  users.  u 

possible  to  obtain  sets  of  coupled  equations  f  ■ r  t  •  •  '  . 

state  equilibrium  distribution  of  these  slates.  T  n  •  »  •  • 

solved  numerically  and  curves  for  the  average  w;, .  1  .  • .  g  , 
service  and  total  delay,  are  obtained.  The  model  ;r,  t  in 
paper  leads  to  one-dimensional  Markov  chi  ins,  as  opposed  to 
two-dimensional  Markov  chains  obtained  in  (1). 

1 1  THE  MODEL 

The  user  terminal  is  shown  in  Fig.1.  The  buffer 


-2- 


15) 


has  unlimited  capacity;  time  is  slotted.  Packets  are 
assumed  to  be  of  fixed  length  L,  with  one  time  slot 
equivilant  to  a  packet  length.  It  is  assumed  si  so  that  the 
sensing  of  a  collision  is  instantaneous.*  The  arrival 
process  is  modeled  as  a  Bernoulli  process  with  average  rate 
of  arrival  (or,  equivilantly  probability  of  an  arrival  in  a 
given  slot)  equal  to  o-  . 

The  packet  of  the  head-of-line  (HOL)  in  the  buffer 
enters  the  transmitter,  with  probability  f,  and  is 
transmitted.  At  the  same  time  a  replica  of  the  packet 
remains  in  the  buffer  (as  HOL).  In  the  case  of  success, 
(the  terminal  is  in  the  unblocked  state),  the  replica  packet 
is  discarded  and  the  following  packet  becomes  HOL.  In  the 
case  of  collision,  (the  terminal  is  in  the  b.loc_ked_  state), 
the  replica  packet  enters  the  transmitter  with  probability 
f,  (the  same  probability  by  which  the  HOL  packet  in  the 
unblocked  terminal  is  transmitted),  for  retransmission  .  A 
packet  that  arrives  at  an  idle  user  goes  directly  to  the 
transmitter  without  incrementing  the  buffer  content. 

Hence,  the  terminal  (or  user)  can  be  in  one  of  two 
states,  an  idle  state  or  an  state.  In  the  idle 
state,  the  terminal  is  unblocked  and  the  buffer  is  empty, 
hence  a  packet  is  generated  with  probability  o  and 
transmitted  instantaneously.  In  the  .votive  state,  the 


*  Instantaneous  sensing  of  collisions  is  clearly  not 
possible  for  satellite  channels. 
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packet  (whether  it  is  the  HOL  packet  in  an  unblocked  state 
or  it  is  the  replica  packet  in  the  blocked  state),  is 
transmitted  with  probability  f. 

In  the  following  we  describe  the  user  Markov  chain 
and  the  system  Markov  chain: 
a)  The  User  Markov  Chain. 

The  state  of  the  user  is  denoted  by  (i),  where 
i  =  0  ,  1  ,  2  ,  .  .  .  indicates  the  number  of  packets  in  the  buffer. 
Let 

7 T.  4  steady  state  probability  of  the  state  (i) 

and 

oo  i 

G(z)  4  v  7Tj  /  (  1  ) 

i=o 

In  order  to  characterize  the  transition 
probabilities  of  the  user  Markov  chain,  we  neea  to  define 
the  following  additional  quantities; 

r  =  Prob. [successful  transmission/user  is  active] 

=  Prob.  [success/user  is  idle] 

=  O'  q 

where  q  is  the  probability  of  success,  assuming  a  packet  has 
already  arrived,  given  the  user  is  idle.  Now,  the 
transition  probability  of  the  user  Markov  chain  can  be 
obtained  as  shown  in  Fig. 2,  whore: 

A  =  Prob.  [no  packets  arrive  to  the  active  terminal  and  a 
packet  succeeds] 

=  (  1  -cr  )  r 

B  =  Prob.  [no  packets  arrive  and  no  transmission,  or 
transmission  and  no  success,  OR  1  packet  arrives  and  one 
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] 

success  ]  i 

=  (  1  -cr  )  (  1  -r  )  +c r  r 

C  =  Prob.  [one  arrival  and  no  transmission,  or  transmission 
with  failure] 

=  cr  (  1  -r  ) 

D  =  Prob.  [1  arrives  and  collides] 

=  cr  (  1  -q  ) 

E  =  Prob.  [no  arrivals  or  1  arrives  and  succeeds] 

=  (  1  -o'  )  +  cr  q 

Hence 

=  A  Vi  +  JiTrj  +C7rj-i  j>2  (2,a) 

n1  =  A rr2  +  l.5jr  +  Dn0  (2.  b) 

n0  =  +  )•:  (  2  .  c  ) 

Applying  Eqn.(l),  we  get: 

G(z)  =  [Az  1 +B  +  Cz  ]G ( z )  +  [ -Az  ^  +  ( D-C ) z  +  (E-B  )  ]  7 T  (  ■*  ) 

o 

Differentiating  (3)  with  respect  to  z  and  evaluating  at  z=1, 

we  get 


7T 


o 


_A-C  _ 
A-C  +  D 


i  .  e  . 


7T  = 

°  r -cr  (  1  -q  ) 

To  determine  the  average  buffer 

differentiate  (3)  twice  and  evaluate  at  z=1; 
✓ 

Q  =  G ( 1 )  as; 


(4) 

size,  0  ,  we 

we  then  obtain 


A  [  1  -  7ro  j 


1  -o'  )rf  1-7T  ] 
_ _ _ o 

r-ir 


b  )  The  System  M nrk o  v  _C h_'ijn  ._ 

The  state  of  the  system  is  described  by 
of  active  users.  Let 

M  =  total  number  of  users  in  the  system 

and 


n  =  number  of  active  users 


Then 


M-n  =  number  of  idle  user.s 

P  ^  4  steady  stat  probability  of  havin 

users  in  the  system. 

Fig.  3  shows  the  system  Markov  chain, 
the  transition  probabilities  we  need  the 
definition; 

g  ( n  )  =  Prob  .  t  i  out  of  n  blocked  use 
transmission ] 

=  ( J)  f 1 ( 1 -f ) n_i 


q  .  (  n  )  = 
J 


P  r  .)  b  .  [  j  out  of  M-n  idle  user 
transmission  ] 


.("TV  j(.-cr  I"-"'-’ 


E  ^  The  conditional  probability  that 


( ‘3  ) 

the  number 


g  n  active 

To  obtain 
following 

rs  attempt 


attempt 


t  li  e  buffer 


size  equals  unity  given  that  the  user  is 


Referring  to  Fig. 3,  we  can  write  P  as; 

>n  =  c|Q(n+  i)^  (  n+1)  E  Pn+i  +  |fi0(n)|l  -^(njK)  +qi(n)gQ(n)J  i'n 


1  <n <(  M -  i  ) 


(6.  a  ) 


=q0(l)Kj(2)  K  I*2  +  {qo(l  )|i-S1(l)K|  +  qi(l)fi0(l)i»>1  <  6  .  b  ) 


=  qo(  1)  gl  (1  )  E  Px  +  |  qQ  (  o)  +qt(  M)  j  PQ 


(  6  .  c  ) 


The  first  term  in  Eqn.  (6. a)  is  the  probability  that 


only  the  active  user,  with  buffer  contents  equal  to  unity, 
transmits,  times  the  probability  the  system  is  in  state 
(n+1).  The  second  term  is  the  probability  that  either  none 
or  at  least  two  active  users  or  one  active  with  buffer 
contents  greater  than  one  transmit,  or-  no  active  but  one 
idle  transmit,  times  the  probability  that  the  system  is  in 
state  (n).  The  third  term  is  the  probability  that  one  idle 
user  and  at  least  one  active  user  transmits,  times  t  hi  e 
probability  of  being  in  state  (n-1).  The  last  term  is  the 
probability  that  the  i  idle  users  generate  and  transmit  a 
packet  times  the  probability  of  being  in  state  ( n  -  i  )  .  Vi  o 
then  sum  over  i , ( 2<i<n ) . 


III.  SUCCK3S  AND  COLLISION  PHOIUfUL  I  T  I KS  . 
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The  success  and  col'  is  ion  probabilities  are 


determined  as  follows; 

r  =  Prob.  [successful  transmission/user  is  active] 

f  V(l-o')M-,1(l-f)n~1  P 

i)  »  1 _ _ 

1  "  Po  ( 7  ) 

Where  the  expression  in  the  denominator  represents  the 
probability  that  there  is  at  least  one  active  user.  The 
numerator  represents  the  probability  that  the  active  user 
under  consideration  attempts  transmission  times  the 
probability  that  (M-n)  idle  users  don't  receive  a  packet 
arrival  and  (n-1)  remaining  active  users  don't  attempt 
transmission  times  the  probability  that  there  are  n  active 
users  in  the  system.  We  sum  over  all  possible  values  of  n. 
Similarly, 

H"1  ,  .  iM-n-l  ,  n 

o'  L  (  1  -v)  ( 1  -  f  )  l> 

(ij  =  n_o  u 


1  -  J> 


M 


=  trq 


(8) 


1 V .  AVERAGE  PACKET  DELAY 

The  average  total  delay  per  packet  is  the  sum  of  the 
average  waiting  time  in  the  buffer,  VI,  plus  the  average  time 
spent  from  the  first  attempted  transmission  to  the  final 
successful  one,  (the  service  time  S). 

From  E q n . ( 5  )  and  applying  Little's  result,  we  get 

W  =  5 

<r 
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=  r  (  1  -q  )_ 

r-cr  q  (  9  ) 

To  determine  the  service  time  .9  ,  wo  note  that  the 


average  transmission  time  given  that  the  terminal  is  either 
idle  and  collision  occurs  or  the  terminal  is  active,  equals 


to 


oo 

S' 

n-  l 


n  (1  -  r) 


n-i 


=  2 

r 


When  success  occurs  on  the  first  attempt  of  transmission, 
the  service  time  is  one  packet  long.  Hence  S  is  given  by 


S  =  —  [  (  1  -  7T  )  +  7T  (  1  -  q  )  ]  +  1  x  7Tq 
r  o  o  o 

=  ft1-7Toq]  +  7Toq 


(10) 


V.  NUMERICAL  results 

The  set  of  coupled  non-linear  equation  (2)  through 
(8)  were  solved  numerically  to  obtain  the  probabilities 
[  R  j  J  t  77^  ,  7^  t  and  Q.  This  has  been  carried  out 
iteratively  as  follows:  Referring  to  Fig.1),  we  start  with 
an  initial  value  for  E  ,  the  conditional  probability  that  the 
buffer  size  equals  unity  given  that  the  user  is  active, 
solve  the  set  of  linear  simultaneous  equations  for  the 
system  Markov  chain  and  obtain  [ P  ] .  With  these  values  of 
Pj’s,  we  obtain  the  values  of  r  and  q  as  given  by  equations 
(7)  and  (8)  and  then  using  Eqns.  ( 2 .  c  )  ,  (9)  and  a  proper 
iteration  scheme,  we  obtain  the  new  value  for  E.  In  Fig. 5, 
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and  for  M  =  4  ,  P  =  0.3,  wo  plot  1 1  <  p r  o b  a b  i  L  i  t  i  o  s  P  P  and  P  ^ 
versus  throughput.  Fig. 6  is  a  plot  for  the  probabilities 
7T0,  7 T.|  ,  and  E  versus  throughput.  Then  Fig.  7  is  the  delay 

throughput  curves.  We  notice  that  for  M  =  U  and  P  =  0 . 3 ,  the 
maximum  throughput  is  0.466,  while  for  the  model  presented 
in  (1)  and  for  the  same  value  of  M  and  P ,  it  is  0.434. 

V 1 .  CONCLUSION 

This  paper  considered  the  problem  of  interacting 
buffered  terminals  in  multiple  access  sc  h  ernes.  The 
analysis  of  the  slotted-ALOH A  multiple  access  scheme  with  a 
finite  number  of  buffered  terminals  lias  been  presented. 
The  analysis  is  based  on  the  idea  of  using  two  Markov  chains 
(system  end  user  Markov  chains)  that  partially  accounts  for 
the  interaction  between  the  users.  The  model  here  is 
simple  and  leads  to  one-dimensional  Markov  chains. 
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Worst  Case  Aquisit'nn  of  PN  Sequences 
Farhad  Ilemmati  Donald  L.  Schilling 
Department  of  Electrical  Engineering 
City  College  of  New  York 
New  York,  N.  Y.  10031 

Abstract 

An  upper  bound  on  the  worst  case  autocorrelation  of  PN  sequences  of  a 
given  length  is  obtained.  The  value  of  the  worst  case  autocorrelation  curve 
at  some  specific  points  is  found  and  using  interpolation  methods  a  smooth 
curve  fitting  to  the  worst  case  autocorrelation  curve  is  formulated. 


Introduction 

Letg(0),  g(l),  g(2),’**>  g(L-l)  be  a  PN  sequence  of  length  L  =2n-l. 
The  autocorrelation  of  the  PN  sequence  is  defined  as 


y-i 

V'n=r 

j=o 


+  e  (j+T) 


(i) 


where  7  is  the  partial  correlation  length.  For  T=  0,  i.c.,  for  correlated 
sequences  PQ  (>)  =7  .  Furthermore,  HT(E)  = -1  for  1K  r  ^  L-l.  Hr(y) 
depends  on  tiie  specific  selected  PN  sequence  and  the  phase  T. 


For  a  given  value  of  r  let  us  define 


h(k)=  g(j)  +  ft(j+T)  (la) 

The  shift  and  add  property  of  PN  sequences  implies  that  h(0),  h(l),  h  (2),  •  ■  • , 
h(L-l)  is  again  a  PN  sequence.  Therefore  Hr  0)  depends  on  the  arrangement 
of  the  zeros  and  ones  in  the  PN  sequence  h(k).  Figure  1  shows  a  typical 


The  research  performed  for  this  paper  was  partially  supported  u  .der 
AFOSIi  grant  05431  and  HA  DC  Contract  F10G28-00-C-0095 
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example  of  the  correlation  of  such  a  PN  sequence. 

The  worst  ease  of  UT(V),  for  all  7  ,  is  the  one  which  over  the  set  of 
all  l’N  sequences  of  a  given  length  and  all  possible  phases  1-  L-l,  is 
as  close  to  H  (y)  as  possible.  A  sequence  h (k)  produces  the  worst  case 
autocorrelation  if  it  starts  with  the  longest  possible  run  of  zeros,  which 
is  n-1,  followed  by  the  shortest  possible  run  of  ones  followed  by  the  longest 
possible  run  of  zeros  followed  bv  the  next  shortest  possible  run  of  ones  and 
so  on.  This  sequence  was  invented  by  Davidovici  and  Schilling  [lj and  is 
known  as  a  pscudo-pseudo-noisc,  (PPN),  sequence.  As  an  example  the  PPN 
sequence  of  length  L=31  is  0000100010010010110110111011111.  The  PPN 
sequence  is  very  close  to  the  actual  worst  case  of  the  PN  sequences  for  short 
length  codes.  Unfortunately,  it  was  observed  that  for  long  codes  the  obtained 
upper  bound  on  the  worst  ease  autocorrelation  of  PN  sequences  based  on  the 
structure  of  PPN  sequence  is  not  very  tight.  In  this  paper  we  add  to  the  above 
construction  procedure  the  condition  that  each  n-luplc  along  the  sequence  must 
appear  only  one  time  (window  property).  This  condition  results  in  the  following 
well  known  construction  technique: 

1.  The  Worst  Case  .Sequence  Generation  Algorithm 
t.  Generate  n  ones  and  set  i  =  0. 

2.  Append  a  zero  and  check  if  the  generated  n-luplc  had  not  previously  occurred. 
If  not  set  i  =  i+1,  go  to  2,  otherwise  change  the  last  generated  bit  to  one,  set 
i  =i  +  l  and  go  to  3. 

3.  Lf  i=2U  go  to  4.  Otherwise  go  to  2. 

4.  Delete  the  first  n+1  bits  and  stop. 

1 

Example:  The  generated  sequence  for  n=f>  is 

niiio  ooooi  ooui  l  ooi  01001 1 1 0101 1 01 1 1 1 1 . 
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The  above  construction  technique  was  originally  given  by  Ford  12J  and  is 
generalized  by  Colo  mb  [ijj  and  is  extensively  studied  by  Fredrieksoi^J- 
Note  that  the  generated  sequence  satisfies  the  run  property  and  the  window 
property  of  the  PN  sequences  while  the  I’PN  se<iuencc  only  satisfies  the  run 
property.  Moreover,  the  sequence  constructed  by  the  alx>ve  algorithm  can 
be  generated  by  an  n-slage  nonlinear  feedback  shift  register. 

A  pure  cycling  shift  register  is  an  n-slage  shift  register  with  the  last 
stage  connected  to  the  first  stage  of  the  register.  The  length  of  the  cycles 
generated  by  the  pure  cycling  shift  register  are  divisors  of  n  and  the  number 
of  cycles  is  [YjJ, 

Z(n)=^I  0(d)2n/d,  (2) 

n  d|n 

where 0  (d)  is  the  Kuler's  0-function.  hi  the  following  we  summarize  the 
properties  of  the  nonlinear  sequences  generated  by  the  above  algorithm  [6j: 

Property  1  ; 

The  nonlinear  sequence  is  constructed  by  joining  the  cycles  of  the  pure 
cycling  feedback  shift  register. 

Property  2 : 

The  breadth  of  a  vector  is  defined  as  the  longest  number  of  contiguous 

zeros  in  the  vector.  Then:  the'  algorithm  exhausts,  in  order,  the  cycles  of 
breadth  n-1,  n-2,  n-3,  etc. 

2.  Formulation  of  the  Upper  Hound 

Properties  1  and  2  staled  ab  vc  can  be  used  to  find  the  value  of  the  worst 
case  autocorrelation  at  some  specific  points.  The  first  exhausted  cycle  is  the 
cycle  of  breadth  n-1.  Hence 

Hr  (n)  =  n-2  (3) 
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There  is  one  cycle  of  breadth  n-2,  two  cycles  of  breadth  n-3,  and  in 
general  2°  ^  cycles  of  breadth  b  for  -  b  -n-2,  where  for 

nockl  and  N  = -^  for  n  even.  The  weight  of  a  cycle  is  the  number  of  ones 
in  ttiat  cycle.  For  I)  there  arc  less  than  2  cycles  of  breadth  b 

since  all  of  the  low  weight  cycles  of  length  n  have  already  been  exhausted. 

Since  a  cycle  of  length  less  than  n  cannot  have  a  breadth  greater  than  Nj, 
,all  the  cycles  exhausted  up  to  the  breadth  arc  of  length  n. 

The  partial  correlation  of  two  sequences  x  and  y  can  lie  also  defined  by 


Rt  (Y)  =  y-  2-Ufo7  (X+  y) 


(4) 


y 

where  W  (z)  is  the  Hamming  weight  of  the  sequence  y  =  x  +  y  truncated  to 
.length  7.  Using  property  2  and  Eq.  (4)  we  can  find  Nj  points  of  the  worst  case 
autocorrelation  by  the  recursive  relation 


Ht  (n  •  2j)  =  UT  (n-21"1)  +  n.2M  -  2  [■’ 


for  1  —  i  -  N 


(•r>) 


The  second  term  is  the  right  hand  sifle  (rhs)  of  (5)  is  the  number  of  bits  in 
all  cycles  of  breadth  b=n-l-i.  The  third  term  in  the  rhs  of  (f>)  corresponds 
to  the  fact  that  in  the  set  of  all  binary  (i- 3)_Luple.s  the  number  of  ones  is 
equal  to  the  number  of  zeros,  and  the  last  term  of  (5)  considers  the  two 
ones  at  the  beginning  and  at  the  end  of  each  breadth. 

Equations  (3)  and  (5)  can  be  rewritten  as 


'r^n)  =  n-2 

^  HT(n*  2*)  =ryn.2  1  *)  t  (n-i-3)25  1  for  1-i-  Nj 


or 


Ht  (n-2*) 


n-2  +  V  (n-j-3)  2' 

J=1 


J-l 


1-  i  — 


(6) 


for 
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*  i  1  j 

Using  v  j«lr  =  (i— 1)2  +1  in  (G),  the  final  result  is 

j=l 

f  Ht  (n)  =  n-2 

[uT(n.2,)  =  (n-i-2)2l  ,  l^i^Nj 


(?) 


Observe  that  in  the  set  of  all  cycles  of  a  given  breadth  b,  the  last 
cycle  exhausted  is  the  cycle  consisting  of  b  zeros  followed  by  n-b  ones. 
Application  of  this  observation  in  (7)  gives  another  set  of  points  of  H  (^j, 
namely 


'  H  (n-1)  =  n-1 

- 

I?r  (n-  21-l-i)  =  (n-i-2)2i+i+l  ,  1-  i  -  Nj  (8) 

nj 

The  formulation  of  RT(7)  for  7  >  n*  2  in  general  is  rather  compli¬ 
cated.  However  for  the  case  that  n  is  a  prime  number,  we  note  from 
(1)  that  the  cycles  of  the  pure  cycling  register  arc  of  length  one  or  n. 

In  this  case  there  arc  two  cycles  of  length  one  and  there  arc  K  cycles 
of  length  n  and  weight  W  where 


1  ,n 


n  'W 


0<  w  -  n- 1 

f 


(9) 


In  order  to  find  an  upper  bound  on  the  value  of  the  worst  case  auto¬ 
correlation  we  assume  that  n  is  a  prime  and  the  worst  case  sequence 
begins  with  the  cycle  of  length  n  and  weight  one,  followed  by  the  cycles 
of  weight  two,  followed  by  the  cycles  of  weight  three,  etc.  With  this 
assumption  the  worst  case  autocorrelation  at  partial  length 7  is  upper 
bounded  by 
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«T  (7) 


—  ^  n*N..n-  2  £  WN" 


w 


w 


W-l 


w 


where  k  is  the  weight  of  the  highest  weight  cycles  in  the  partial  length 
^  .  Substituting  (9)  into  (10)  the  result  is 


Rr  (7) 


k  k 


£}\>?~2<L  n  \v' 

W-l  W- 1 

W  ,n  .  /  n-1  \  . 

Since  ~(w)  =(\v_  1 1  we  have 

k  k 

O'2-  Sw’i' 

W-l  W  \v-lw_1 


(10) 


(11) 


where 


k  k 

T  =  n  I'  Nn  =  T  (ni 
w  *  ->  'w 
W-l  W=1 


(12, 


In  particular  for  k  =  (n-1)  wc  have  from  (12) 

n-1 


2 

7  =  v 
W=1 


0 


and  from  (11) 

n-1 

2 

nT  (2n-1-i)*  I' 

W-l 


c 


n-1 
2 

)  -  2  v 

w-i 


C>  -  c‘> 


(13) 


Using  Stirling's  Formula  for  factorials,  (13)  reduces  to 


HT(2n-J_l)<  2 11 

-17 

For  n  large,  1/  ~  2n  -  2  w  2°  and 

I{T  (L'/2)  -  - ~ - 

\S'z  W(n  -  1) 


s 


(14) 
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To  find  a  smooth  curve  for  the  upper  boim<  1  on  the  worst  case 
autocorrelation  of  the  PN  sequences  we  can  use  any  interpolation 
technique,  c.g.  I-ngrange's  method  among  the  specific  values  of 
Kt  (Y)  given  in  (7),  (8),  (11),  and  (14).  The  simplest  case  is  an 
interpolation  among  the  three  points  (0)  =  0,  ltT  (L')  =  0  and 
H_(L'/2)  given  by  (14),  which  yields  a  second  degree  polynomial, 

Rt  (Y)  **  - - (I-y/L*  )  (15) 

i 2 it  (n-1) 

Obviously  a  higher  degree  polynomial  will  give  a  better  approxima¬ 
tion  to  the  worst  case  autocorrelation  curve. 

hi  Fig.  2  Ihc  worst  case  curve  for  n=ll  generated  by  the  algorithm 
described  in  Sec.  1,  is  compared  with  the  approximation  curve  of  (15) 
and  a  five  point  interpolation  curve.  Figure  3  is  the  Same  as  Fig.  2  but 
for  n=13.  It  should  be  mentioned  that  in  the  formulation  of  HT(  -,)  we  have 
assumed  that  the  worst  case  curve  is  symmetrical  with  respect  to  the  line 
Y  =  L'/2. 

In  the  case  that  n  is  not  a  prime  number  the  length  of  the  cycles  of 
the  pure  cycling  shift  register  arc  divisors  of  n.  An  upper  bound  on  the 
worst  case  autocorrelation  curve  for  n  a  nonprime  number  can  again  be 

obtained  similar  to  the  upper  Ixmnds  or  (11)  and  (12).  Since  the  average  weight 

per  length  of  the  short  length  cycles  is  approximately  equal  to  the  average  weight  per 

length  of  the  long  length  cycles,  the  upper  bounds  of  (11),  (12),  and  (15) 
obtained  for  n  a  prime  number  are  valid  for  any  n. 

Figures  4  and  5  compare  the  worst  case  autocorrelation  curve 
generated  by  the  algorithm  for  n=8  and  n=10  with  the  three  point  interpolation 
curve  of  Fq,  (15). 
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3 .  Conclusions 

This  paper  has  presented  an  upper  bound  on  the  worst 
case  partial  autocorrelation  of  nonlinear  PN  sequence. 

The  results  take  the  form  of  a  readily  implementable 
algorithm  and  an  approximate,  easily  used,  equation  which 
upper  bounds  the  algorithm,  namely  (15). 

It  is  interesting  to  compare  (15)  with  the  results 
obtained  by  Davidovici  and  Schilling  [1J  who  found 

R  ( Y ) <  Y  (1  -  ,  L  =  2n-l  (16) 

T  L 

Note  that  both  results  are  of  the  same  form,  and  partial 
autocorrelation  of  [l ]  exceeds  the  bound  developed  in  this 
paper  by  the  factor 


F  = 


(17) 


Hence  for  n=13,  F  =  1.77. 

When  using  (15)  or  (16)  to  determine  the  probability 
of  false  alarm  or  detection  when  acquiring  a  PN  sequence, 
a  value  of  F  of  the  order  of  (17)  is  often  extremely 
significant. 
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.  4  Worst  case  si-puenco  of  length  l.~  2.r>ri  nencrnU'd  l>v  ( lu*  algorithm  anti 
a  ( h ift*  point  intorpolnlion  curve. 


Fig.  5  Worst  case  sequence  of  length  L  -  1023  generated  by  the  algorithm  and 
interpolation  curve. 
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Summary 

An  infinite  cyclic  sequence  of  period  p,  where  p  <,2 

for  some  positive  integer  k,  is  a  shift  register  sequence  if 

any  block  of  k  consecutive  bits  over  one  period  of  this 

sequence  is  not  repeated.  The  sequence  of  period  p  =  2  -1  is 

known  as  the  maximal  length  sequence  and  the  sequence  of 

length  p  =  2k  is  called  a  deBruijn  sequence  [1].  For  a  given 

2k“1-k 

k  the  number  of  deBruijn  sequences  is  2 

A  linear  maximal  length  sequence  (also  known  as  a  PN 
sequence)  can  be  generated  by  a  k-stage  linear  feedback 
shift  register  (LFSR)  with  a  primitive  characteristic 
polynomial  [2].  For  a  given  k  the  number  of  linear  maximal 
length  sequences  is  0(2  -1)/k  , where  0  is  the  Euler's 
0-function.  An  efficient  method  to  construct  the  nonlinear 
sequences  is  the  cycle  joining  method  [2], [31  which  is  based 
on  joining  the  sequences  of  short  period  together  to  obtain 
a  deBruijn  sequence. 

In  the  following  a  new  property  of  the  PN  sequences  is 
presented.  For  generation  of  nonlinear  sequences  of  period 
p  =  2  this  property  and  several  other  well  known  properties, 
[2],  of  the  PN  sequences  are  used  to  join  the  cycles  of  an 
LFSR  with  characteristic  polynomial  G  (  x  )  =  g  (  x  )  ( 1 +x  ) 11 ,  where 
g(x)  is  a  primitive  polynomial  of  degree  m  and  G(x)  is  of 
degree  k=m+n. 

A  block  of  length  m+j  consecutive  bits  on  a  PN  sequence 
generated  by  an  m-stage  LFSR  is  called  an  extended  state  of 


degree  j  of  the  shift  register  and  is  represented  by  the 
(m+j ) -tuple  S=(sm+j ,sm+j_1 , . . . , s2,s 1 ) .  The  conjugate  of  S, 

C  3  ].,  is  defined  by  S  =  (  sm  +  .  ,  sm  +  .  _  1 . s^s^,  where  s1  is 

the  binary  complement  of  s^.  A  characteristic  state  over  a 
PN  sequence  is  defined  as  an  extended  state  of  degree  j  over 
the  PN  sequence  such  that  its  conjugate  satisfies  the 
recurrence  relation  of  (1+x)1-*. 

Property :  There  exists  a  characteristic  state  on  every  PN 
sequence  for  any  j  s  1  . 

Proof :  The  proof  is  by  induction.  The  property  is  obviously 
true  for  j  =  1.  Since  there  exists  a  block  of  m  consecutive 
ones  followed  by  a  zero  ,  (  1  ,  1 ,  .  .  .  ,  1 , 0 )  ,  on  every  PN  sequence 
and  the  conjugate  of  this  extended  state  is  a  block  of  m+1 
consecutive  ones  (1,1,...  ,1,1)  which  satisfies  the 
recurrence  relation  of  1+x.  Now,  if  we  assume  that  the 
property  is  true  for  j=i,  Lempel's  homomorphism  [3]t  can  be 
used  to  prove  that  the  property  is  also  true  for  j  =  i  + 1 . 
Moreover,  the  proof  gives  an  itterative  method  to  find  the 
characteristic  states  of  a  PN  sequence  for  any  j  a  1  . 
Characteristic  states  of  the  PN  sequence  1110010  for  1<j<8 
is  shown  in  Table  I. 

Let  be  the  set  of  all  cycles  of  a  LFSR  with  the 
characteristic  polynom*  '  g( x ) ( 1 +x ) n , [ 4 ] .  Ck  can  be  divided 
into  four  subsets  A,B,C  and  D.  The  subset  A  consists  of 
cycles  of  length  (2m-1)2e^^  which  satisfy  the  recurrence 


relation  of  g(x)(1+x)J,  0^  jin-1,  where  e(j)  is  an  integer 
such  that  2e ^ ^ ~ 1  < j <  2e ^ ^  .  The  subset  B  consists  of  cycles 
of  length  ( 2m- 1 ) 2e ^ n ^  which  satisfy  the  recurrence  relation 
of  g(x)(1+x)n.  Similarly  the  subset  C  consists  of  cycles  of 
length  2e  ^  which  satisfy  the  recurrence  relation  of  ( 1 +x  )  ^ 
0^j^n-1.  And  finally  the  subset  D  consists  of  cycles  of 
length  2ev  '  which  satisfy  the  recurrence  relation  of 
(1+x)n.  The  existence  of  a  characteristic  state  on  every  PN 
sequence  for  any  j  ^  1  can  be  used  to  prove  the  following 
theorem . 

Theorem :  Let  a  be  a  cycle  of  length  (2m-1)2e^  in  the  set 
A.  There  exists  2e ^ ^  states  on  a  such  that  their 
conjugates  are  in  the  set  D.  The  conjugates  of  the 

( 2m-2 ) 2e ^ ^  remaining  states  on  a  are  in  the  set  B. 

Corollary :  Let  0  be  a  cycle  of  length  (2m-1)2e^n^  in  the 

set  B.  There  exists  2e  states  on  /S  such  that  their 

conjugates  are  in  the  set  C.  The  conjugates  of  the 

(2m-2)2e^n^  remaining  states  on  fi  are  in  the  set  A. 

Since  the  number  of  cycles  of  an  LFSR  with  the 
characteristic  polynomial  g(x)(1+x)n  increases  exponentially- 
with  respect  to  n  ,[3],[4],  the  number  of  different  deBruijn 
cycles  which  can  be  constructed  by  joining  the  cycles  of 
g(x)(1+x)n  increases  double  exponentially  with  respect  to  n. 
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Table  I.  Characteristic  States  of  the  PN  Sequence 


1110010  for  1  <  j  <  8. 


s 

A 

s 

( 1+x)^ 

1110 

1111 

1  +x 

0101  1 

01010 

1.x2 

1  10010 

110011 

1+x+x2+x3 

1011100 

1011101 

1.x4 

100101 1 1 

10010110 

4  S 

1+x+x  +x 

01 1100101 

01 1100100 

.  2  4 

1+x  +x  +x 

00101 11001 

00101 11000 

2  1 

1+x+x  +x 

1 1 100101 110 

11100101111 

1.x8 
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Superresolving  image  restoration  using  linear  programming 


R.  Mammone  and  G.  Eichmann 


Snpcrre-»olviiij:  image  restoration  tSIU)  in  the  preseme  «>|  noise  i'  « »*iisnlere<l  IVw  SIR  alj'-n  it  Inns  h.i\i- 
•  leu  must  rated  the  abiht  v  to  lesolve  two  point  soiiri  es  sp.n  ed  one  lull  ol  I  In-  R.t\  Iri;;  h  dist  a  nee  ap.n  l .  In  this 
paper,  it  is  shown  that  the  SIR  <>l  a  two  point  noncoherent  sou  it  *•  spaced  one  tenth  « »f  a  Wavleigh  distance 
apart  is  possible.  'The  met  hod  presented  uses  opt  i  mat  data  lining  leclinitpies  lia^t-d  on  the  inelliods  ol  linear 
programming.  For  noisy  images,  a  combination  of  linear  eigenvalue  prefiltering  and  optimal  data  lilting 
is  usi*il .  It  is  also  shown  that  (nr  a  diffract  ion-limited  image  of  two -point  sources  spaced  one  halt  ol  the  Ray¬ 
leigh  distance  apart,  where  the  input  is  contaminated  with  significant  noise,  SIR  is  achievable.  'These  re¬ 
sults  have  important  implication^  m  atmospheric  physics,  geoplp  sics.  ia«iio  aslioiiomv,  medical  diagnostics, 
ami  digital  bandwidth  compression  applications  where  the  decnnvohiliou  of  noisy  bandwidth  compressed 
images  is  one  ot  t  tie  fund  a  mental  limit  at  ions.  The  lechuiipies  descrihefl  ale  speeil  leally  designed  lor  i  in  pul 
sive-tvpe  images. 


I.  Introduction 

"The  problem  of  restoring  a  linearly  degraded  image 
has  been  the  subject  of  extensive  investigations.1  The 
mathematical  formulation  of  the  image  restoration 
problem  is  common  to  many  diverse  fields  such  as  at¬ 
mospheric  physics,-  geophysics,11  spectroscopy,’  radio 
astronomy, as  well  as  the  general  subject  of  numerical 
analysis."  In  general,  the  inversion  of  a  linear  degra¬ 
dation  can  he  accomplished  in  a  number  of  wavs.1 
When  the  measurement  error,  which  is  the  noise,  is 
negligible  and  the  degradation  process  is  well-behaved, 
the  direct  inversion  of  the  degradation  process  can  easily 
he  performed.  For  a  well  conditioned  degradation 
process,  in  the  presence  of  noise,  methods  such  as 
Wiener  filtering  are  appropriate.11  However,  in  many 
cases  the  degradation  process  is  ill-conditioned.  In  this 
case  the  measurement  errors  are  greatly  amplified  in  the 
restoration  process  leading  to  an  image  estimate  that 
is  dominated  by  noise.  In  this  paper  we  discuss  the 
restoration  of  ill-conditioned  linearly  degraded  images 
that  have  been  corrupted  by  appreciable  noise. 

While  many  types  of  ill-conditioned  image  restoration 
can  he  considered,  a  particular  ill-conditioned  image 
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restoration,  the  restoration  of  a  diffraction  limited 
image  of  finite  extent,  has  been  used  in  this  paper  as  a 
bench  mark  to  provide  a  stringent  test  of  the  algorithm. 
Furthermore,  this  model  has  many  other  important 
applications  such  as  the  restoration  of  a  severely 
bandwidth-compressed  video  images.  Such  images  are 
of  interest  in  digital  TV  applications. 

It  has  been  shown'-'-1'1  that,  for  images  of  finite  extent, 
in  the  absence  of  noise,  spectral  components  of  the 
signal  that  have  been  removed  can  he  reconstructed 
from  the  low  passband  information  using  analytic 
continuation.  The  extrapolation  in  the  spatial  fre¬ 
quency  domain  represents  an  increase  in  the  spatial 
resolution  in  the  spatial  domain.  An  image  restorat  ion 
process  which  provides  this  increase  in  resolution  is 
called  a  superresolving  image  restoration  (SIRl  process. 
A  complete  survey  of  recent  restoration  methods  in¬ 
cluding  SIR  has  been  given  by  Frieden.1  In  general, 
statistical  image  restoration  methods  are  not  superre¬ 
solving.  Deterministic  methods  use  various  forms  of 
regularization11  to  provide  stability  for  the  image  res¬ 
toration.  Linear  shift-invariant  filters,  while  they 
provide  stability,  are  also  not  superresolving.  It  lias 
been  postulated  that  only  nonlinear  filtering  methods 
will  provide  SIR  in  the  presence  of  noise.1  The  purpose 
of  this  paper  is  to  show  that  the  combination  of  linear 
shift-invariant  filtering  of  the  image  and  nonlinear 
restoration  does  provide  SIR  capability  in  a  realistic 
noise  environment. 

The  optimal  data  fitting  step  uses  linear  programm¬ 
ing  (LPl  techniques.  The  LP  method  provides  both 
cost  and  time  effective  ways  to  produce  SIR.  Although 
methods  similar  to  LP  have  been  suggested  for  image 
processing,1-  11  we  have  specifically  addressed  ill-con- 
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ditioned  problems  and  have  incorporated  a  stabilizing 
presmoothing  step,  which  also  yields  computational 
advantages. 

II.  Formulation 

The  SIR  problem  requires  the  inversion  of  a  Fred¬ 
holm  integral  equation  of  the  first  kind  with  a  sine 
function  kernel.  The  integral  equation  corresponding 
to  the  coherent  1-1)  case  is 


where  a' U  )  is  the  diffract  ion-limited  (Did  image,  /(v) 
is  the  original  image  of  length  2b,  and  /,  is  the  cutoff 
frequency  of  the  transfer  function.  It  is  well-known 
that  the  incoherent  handlimiting  kernel  has  twice  the 
spatial  cutoff  frequency  as  compared  with  the  corre¬ 
sponding  coherent  kernel.  Therefore,  the  coherent  SIR 
is  more  stringent  and  thus  provides  a  better  bench  mark 
than  the  corresponding  incoherent  kernel. 

The  discretized  form  of  the  image  restoration  prob¬ 
lem  can  he  formulated  as  a  solution  to  the  matrix  linear 
equation 

g  =  IK  +  n,  (2) 

where  f.g  are  the  sampled  values  of  the  restored  (de¬ 
sired)  diffraction-limited  and  measured  handlimited 
images,  respectively,  n  is  the  sampled  noise  vector,  and 
//  is  a  matrix  representing  the  sine  integral  operator 
obtained  by  using  some  appropriate  quadrature  rule. 
Without  loss  of  generality  we  shall  take  f,  g,  and  n  to  be 
/V-dimensional  vectors.  Thus  //  is  an  N  X  N  square 
matrix.  The  system  Of  linear  equations,  represented 
by  Kq.  (2),  is  underdetermined  since  there  are  At  equa¬ 
tions  and  2N  unknowns.  The  unknowns  are  the  com¬ 
ponents  of  the  Af -dimensional  vectors  f  and  n.  Ne¬ 
glecting  the  noise  vector  does  not  per  se  simplify  the 
problem  because  for  SIR  the  //,  the  degradation,  matrix 
is  nearly  singular.  Because  of  the  rank  deficiency  of  the 
H  matrix,  the  inverse  H  is  highly  unstable. 

To  stabilize  the  inverse  H  mat  rix,  the  eigenvalues  and 
the  eigenvectors  of  the  inverse  matrix  are  filtered.  The 
eigenvectors  of  the  H  matrix  are  similar  to  the  discrete 
prolate  spheroidal  wave  seq lienees,1  f>  the  analog  of  the 
prolate  spheroidal  wave  functions  of  the  continuous 
case.  The  eigenvalues  of  the  //  matrix  change  abruptly 
from  approximately  unity  to  nearly  zero  as  a  function 
of  the  space  bandwidth  product.  The  high  order  ei- 
gensequenees.  corresponding  to  rapidly  oscillating 
functions,  represent  the  eigenvalues  close  to  zero.  The 
instability  of  the  inverse  matrix  is  due  to  the  near  zero 
eigenvalues  of  the  II  matrix.  By  filtering  (attenuation) 
of  the  image  components  corresponding  to  those  ei¬ 
genvalues,  stable  inversion  can  be  obtained.  There  are 
a  number  of  linear  filtering  schemes  that  have  been 
suggested  to  do  this  filtering. Ir'-17  While  these  schemes 
stabilize  the  inversion  process  of  the  H  matrix,  since 
they  do  so  at  the  expense  of  the  high  order  eigenvalue 
components,  which  are  precisely  those  components  that 
contain  the  high  spatial  frequency  information  neces¬ 
sary  for  supcrresolul ion,  these  methods  are  not  SIR 
methods.  Furthermore,  since  this  approach  does  not 


take  into  account  the  measurement  errors,  a  critical 
factor  in  a  highly  unstable  matrix  inversion,  other 
met  hods  must  be  sought. 

Additional  information  or  constraints  must  be  uti¬ 
lized  toobtain  SIR.  Methods ol  image  restoration  that 
impose  non -negativity  constraints  on  the  restored  image 
to  obtain  stability1-11'  have  been  described.  We  use 
non-negativity,  boundedness,  and  derivative  constraints 
together  with  eigenvalue  filtering  to  obtain  SIR.  The 
addition  of  these  constraints  stabilizes  the  restoration 
process.  This  method  of  restoration  can  be  formulated 
as  a  data  fitting  problem  where  it  is  desired  to  obtain  a 
restored  image,  which,  when  degraded,  is  as  close  as 
possible  in  some  sense  to  the  measured  image.  Two 
measures,  or  norms,  of  closeness  were  invest  igated,  the 
/ 1  and  / ,  norms.  The  / 1  norm  of  a  vector  is  the  sum  of 
t  he  absolute  values  of  the  elements  of  the  vector,  while 
the  /.,  norm  of  a  vector  is  the  maximum  absolute  value 
of  the  vector.  The  method  which  minimizes  this  norm 
is  frequently  called  the  minimax  method. 

Thus  the  method  consists  of  two  steps:  (1)  pres¬ 
moothing  the  measured  image  to  attenuate  or  eliminate 
image  components  which  effectively  could  be  due  only 
to  noise;  (2)  selection  of  the  closest  restored  image  es¬ 
timate  corresponding  to  the  measured  image.  Two 
measures  of  closeness  were  examined  yielding  the 
minimal  l\  and  minimax  methods. 

A.  Presmoothing  the  Measured  Image 

A  method  of  stabilizing  the  problem  is  to  smooth  the 
elements  of  the  measured  image  vector  g.  Ideally  this 
smoothing  can  be  accomplished  by  passing  the  mea¬ 
sured  image  g  through  a  filter  which  eliminates  the 
frequency  components  for  frequencies  higher  than  f, . 
This  filters  out  the  components  due  only  to  noise.  This 
preprocessing  step  is  straightforward  if  the  1)1.  image 
is  provided  in  its  entirety.  However,  the  original  image 
must  be  of  finite  duration  for  a  solution  toexist.  Thus 
the  1)1,  image  must  be  infinite  in  extent,  although,  for 
practical  reasons, some  finite  interval  of  the  1)1,  image 
must  be  selected.  Because  the  measured  image  is  not 
handlimited,  some  other  linear  filtering  process  must 
be  used. 

A  method  of  smoothing  I  he  measured  1)1.  image  can 
be  obtained  by  examining  the  singular  value  decom¬ 
position  (SVl))17  of  the  degradation  matrix  H .  From 
Kq.  (2),  neglecting  the  noise  term,  we  have 

K  = //f  =  ('A,"-Trf.  CU 

where  l  '  and  \’;  are  orthogonal  matrices,  A  is  a  diago¬ 
nal  matrix  comprised  of  the  eigenvalues  of  the  svm- 
mcl  l  ie  matrix  llll 1 ,  and  superscripts  T  and  */.*  denote 
matrix  transpose  and  square  root,  respectively.  The 
standard  form  of  the  SVl )  is  used,  where  the  eigenvalues 
are  in  decreasing  order  from  left  to  right  in  A.  These 
eigenvalues  fall  abruptly  to  practically  zero,  indicating 
I  he  ill-conditioned  nature  of  //.  This  SVl)  decompo¬ 
sition  is  commonly  interpreted  as  a  decomposition  of 
the  image  into  eigenimages.  The  higher  order  eigen¬ 
values  correspond  to  the  higher  order  eigenimages.  The 
amount  of  zero  crossing  is  proportional  to  the  order. 
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Thus  the  high  order  components  correspond  to  the  high 
frequency  content  of  the  image.  We  may  rewrite  Kq. 
(2)  in  a  transformed  coordinate  system 

G  =  A'-F,  (4) 

where 

G  =  U~'g,  (5) 

F  =  Vr  f.  (6) 

The  interpretation  of  G  and  F  is  similar  to  that  of  a 
discrete  Fourier  transform  (DKT).1H  Since  the  ele¬ 
ments  of  the  diagonal  matrix  decrease  abruptly  to  zero, 
it  is  seen  from  Kq.  (4)  that  the  high  order  components 
(high  frequency  components)  are  greatly  attenuated. 
This  is  to  be  expected  since  filtering  in  the  DPT  domain 
is  consistent  with  low  pass  filtering  in  the  diagonalized 
coordinate  system. 

The  smoothing  is  accomplished  by  preprocessing  the 
measured  image  via  the  singular  value  transform  (SVF) 
Kq.  (5)  and  forming 

&  =  U(i.  (7) 

where  G  is  G  modified  by  replacing  all  elements  of  index 
greater  than  r  with  zero.  The  index  r  corresponds  to 
the  largest  order  eigenvalue  of  the  matrix  HHr  whose 
value  is  greater  than  the  standard  deviation  of  the  error. 
This  elimination  of  the  high  order  noisy  component's  of 
the  observed  image  also  yields  a  reduction  in  the  rank 
of  the  degradation  matrix.  'Phis  resulting  rank  defi¬ 
ciency  has  two  effects  on  the  problem:  first,  it  yields 
an  underdetermined  system  of  equations  with  many 
possible  solutions.  We  shall  define  an  optimal  solution 
and  present  a  method  of  obtaining  this  solution  later. 
The  other  effect  is  to  make  the  system  inconsistent;  that, 
is,  the  left-hand  side  of  F,q.  (2),  neglecting  the  noise 
term,  may  have  two  identical  adjacent  rows  of  H  due  to 
the  linear  dependence,  but  as  a  result  of  the  random 
noise,  the  leftTiand  side  will  be  very  different.  Thus 
inconsistencies  in  the  equations  result  from  the  noise. 
Therefore,  we  eliminate  all  but  r  linearly  independent 
equations  in  H.  The  resulting  reduced  system  of  r 
equations  is  well-conditioned  and  consistent.  This 
reduction  of  the  number  of  equations  gives  the  addi¬ 
tional  advantage  of  reducing  the  number  of  computa¬ 
tions  necessary  in  the  next  step. 

At  this  point  the  pseudoinver.se  solution  might 
suggest  itself.  The  problem  with  this  or  any  other  linear 
solution  is  that  it  will  not  extrapolate  the  high  order 
components.  Thus  we  find  a  constrained  optimal  so¬ 
lution  from  this  reduced  system  of  equations  using 
Id*. 

B.  Constrained  Estimate 

'File  general  LI’  problem  can  be  formulated  as  fol¬ 
lows13: 

minimize  cTx,  IS) 

subject  to  Ax  >  b,  (9) 

whrr<*  x  >  0,  <  10) 

and  c  is  called  the  cost  vector,  x  is  I  he  A/ -dimensional 


vector  denoting  the  variables,  the  r  X  M  matrix  A  is 
called  the  constraint  matrix,  and  b  is  an  r -dimensional 
constraint  vector.  The  most  frequently  used  algorithm 
to  solve  linear  programming  problems  is  the  simplex 
algorithm.  This  algorithm  produces  one  of  three  mu¬ 
tually  exclusive  results:  first,  an  indication  that  the 
problem  is  infeasible  (the  problem  is  feasible  if  all  the 
constraints  can  be  simultaneously  satisfied  for  some 
choice  of  variables;  the  problem  is  infeasible  otherwise); 
second,  an  indication  that  the  problem  is  feasible  but 
the  optimal  solution  is  unbounded;  or,  third,  an  optimal 
solution. 

The  minimal  lt  formulation  of  the  image  restoration 
will  now  be  addressed.  We  define  two  /V-dimensional 
vectors  n+  and  n_  so  that  they  contain  the  negative  and 
positive  components  of  n,  respectively.  For  any  given 
row,  one  of  the  two  vectors  n+  and  n~  has  one  zero  ele¬ 
ment,  and  the  other  has  a  positive  element.  Thus  we 
may  write 

B  =  //f  +  n*  -  n",  (11) 

so  that  n  +  ,  n",  and  f  are  all  positive  jV -dimensional 
vectors.  If  we  consider  the  positive  sum  for  row  i, 

n'  +  =  |c,  -  hj\.  (12) 

where  h,  is  the  ith  row  of  //.  This  is  seen  to  be  the  ab¬ 
solute  veVie  of  the  residual  of  the  difference  between  g 
and  bf.  Falling  the  sum  over  all  A’  rows  yields  the  l\ 
norm  ol  t  he  residual  error.  This  will  be  the  cost  func¬ 
tion  that  will  he  minimized  subject  to  the  constraints. 
The  LP  formulation  of  the  minimal  l\  problem  is: 

minimize  /),*  +  ri  J"  +  n*  +  nj  +  .  . .  a*  +  riF  (131 

subject  log  -  //f  )  nf  -  n",  (l.ta) 

wit  h  0  £  f  <  (13b) 

An  alternate  interpretation  of  minimizing  the  resid¬ 
ual  can  be  obtained  by  examining  the  process  in  terms 
of  the  solution  vector.  The  1 1  norm  is  bounded  by  zero. 
If  this  bound  is  obtained,  the  residual  error  is  zero  and 
the  solution  vector  is  completely  in  f.  This  case  per¬ 
tains  to  a  direct  inversion  of  the  degradation  matrix  H. 
in  general,  this  inversion  is  not  possible,  since  the 
measured  image  g  is  usually  not  consistent  with  all  the 
constraints.  Thus  n  will  contain  nonzero  elements,  and 
f  must  contain  corresponding  zero  elements  as  well. 

Following  a  similar  procedure,  the  minimax  estimate 
can  be  obtained  by  LP.  First,  we  let  E  represent  a 
scalar  variable  denoting  the  largest  positive  deviation 
of  the  degraded  estimate  from  the  measured  image. 
Then  we  write  the  LP  formulation  as: 

minimize/'.',  (I  la) 

subject  tn  R  =  lit  +  Ec„  (14b) 

K  =  H  f-Kc,.  (14c) 

where  cr  is  the  r-dimensional  unity  vector.  These 
equations  can  he  solved  liv  the  method  of  LP.  Thus  the 
minimax  formulation  requires  N  more  constraints  than 
the  minimal  /)  formulation  of  the  same  prol  'em. 
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Fig.  1 .  I  >ir**t't  invcrsi*  solul  ion  tor  two  iiu|m]st's  separated  by  one- 
eighth  of  the  Rayleigh  distance. 


Kin.  -  Const  rained  mini  oval  1 1  solution  for  two  impulses  separated 
by  one  eighth  of  the  Kavleigh  distance. 

C.  Numerical  Results 

A  fit'll  re  <d  merit  of  t  ho  SIR  is  t  ho  Kay  hut'll  distance. 
The  Rayleigh  dislitnce  is  inversely  proportional  to  the 
hantlwidl  h  of  t  he  I  muster  fund  ion.  Very  few  restora¬ 
tion  methods  achieve  resolution  of  one-half  of  the 
Rayleigh  distance.1  This  measure  of  resolution  will  be 
used  for  both  noncoherent  and  coherent  i magi nu  sys¬ 
tems.  It  should  he  remembered  that  the  Rayleigh 
distance  lor  a  coherent  imaging  system  is  a  more  strin- 
gont  criterion  than  for  noncoherent  imaging  system. 

The  methods  were  simulated  using  the  simplex  linear 
programming  subroutines  contained  in  the  Interna¬ 
tional  Mathematical  and  Statistical  Library  (IMSId 
KOUTKAN  callable  subroutine  package.  CUNY's 
Central  <  'umpiiting  facility  was  used.  The  bench  mark 
example  is  t  he  dil  l  rad  ion  limitation  of  a  coherent  1  -  I  > 
system,  figure  I  illustrates  reconstruction  of  two  im¬ 
pulses  spaced  one  eighth  of  the  Rayleigh  distance  with 
no  added  measurement  noise.  The  image  is  composed 
of  thirty-two  samples.  The  direct  inverse  solution 
(dotted  line)  demonstrates  very  unstable  behavior,  even 


though  the  only  error  is  that,  due  to  the  finite  precision 
of  the  computer  used  and  arithmetic  errors  which  are 
propagated  numerically  through  the  inversion  process. 
The  exact  solution  is  shown  liv  the  two  impulses  (solid 
line),  figure  2  illustrates  a  perfect  restoration  of  two 
impulses  separated  by  one-eighth  of  the  Rayleigh  dis¬ 
tance  by  the  minimal  l,  norm  method.  The  dotted 
curve  is  the  1)L  image.  It  is  approximately  the  main 
lobe  of  a  sine  function  centered  at  the  center  of  the 
graph.  The  method  completely  restored  the  two  im¬ 
pulses.  The  estimate  is,  the  chain  dot  curve,  completely 
coincident  with  the  exact  curve  (solid  line). 

The  minimax  method  did  not  provide  the  high  reso¬ 
lution  performance  obtained  with  the  minimal  /, 
method  for  two-point  restoration.  Figure  3  gives  an 
indication  of  this  difference.  The  minimax  estimate 
(chain  dot)  curve  is  shifted.  The  second  impulse  is 
down  by  HOT,  and  the  intcrimpulse  trough  is  down  to 
50%  of  the  actual  peak  amplitude.  There  is  also  a  small 
artifact  which  appears  just  to  the  left  of  the  first  im¬ 
pulse.  Although  this  estimate  might  actually  he  ac¬ 
ceptable  for  many  applications,  it  is  clear  that  the 
minimal  /]  estimate  is  preferable.  This  was  found  to 
lie  the  case  in  all  situations  investigated.  The  reason, 
we  believe,  for  this  discrepancy  is  the  increase  in  com¬ 
putations  needed  for  the  ininimax  method  and  thus  the 
greater  the  propagations  of  precision  errors. 

The  method  of  restoration  is  robust  to  noise  yielding 
estimates  with  an  SNR  of  comparable  magnitudes  to 
the  SNR  of  (he  measured  signal.  Here  the  SNR  of  the 
estimate  is  defined  as 

V 

if,  estimate  -  (,  original 

SNH/  -  io  i«»k  — - .  n:»i 

,\ 

L  (/..  ordinal)*' 

i  =  i 

and  the  SNR  of  the  measurement  is  defined  as 


Kijj.  .1.  Constrained  minima*  solution  for  two  impulses  separated 
hv  one  eighth  <»l  tin*  Kayleij;l>  distance. 
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Kin  t.  1 1  estimate  with  no  SVI)  prefilt pring  impulses  separated  l>v 
one-half  of  the  Rayleigh  distance  with  while  noise  added.  <r  =  li.li. 


Ki”.  f>.  1 1  estimate  with  SVI)  pref'iltering,  of  two  impulses  separated 

l>y  one-half  of  the  Rayleigh  distance  with  white  noise  added,  o  = 
<>.:». 

Figure  4  illustrates  the  minimal  / 1  estimate  (chain  dot 
curve)  with  no  SVI)  prefiltering  of  two  impulses  sepa¬ 
rated  by  one-half  of  the  Rayleigh  distance.  White 
Gaussian  noise  with  zero  mean  and  a  n  =  0.2  has  been 
added  to  the  measured  image.  The  corresponding  SNR 
of  the  measured  image  is  21  dB.  The  restored  image 
had  an  SNR  of  14  dB.  The  diffraction-limited  image 
(dotted  line)  has  a  jagged  appearance  due  to  (he  added 
noise.  Figure  f»  illustrates  the  effect  of  the  SVI)  pre¬ 
filtering.  The  example  is  identical  to  that  in  Fig.  4  ex¬ 
cept  that  prefiltering  of  the  measured  image  has  been 
performed.  The  diffraction-limited  image  (dotted  line) 
appears  much  smoother  than  in  Fig.  4.  The  improve¬ 
ment  in  (he  estimate  is  obvious.  The  estimate 
(chain  dot  curve)  lies  almost  completely  coincident,  to 
the  exact  image  (solid  line)  curve.  The  peak  values  of 
the  estimate  are  10%  lower  (han  the  exact  solution,  and 
the  interimpulse  trough  is  10%  higher  than  the  exact 
solution. 

The  method  presented  will  find  r  nonzero  elements 
for  the  2 N  variables.  Thus  the  formulation  given  by 
K<|s.  (Ill)  and  ( 14)  applies  to  impulsive-type  images.  A 


more  general  class  of  images  can  be  restored  if 
smoothness  constraints  of  the  form/,  ()  -  /,  <  of,  are 
appended  to  the  formulations,  where  of,  is  a  bound  on 
the  first  order  difference  (derivative)  of  the  image  to  lie 
restored.  Higher  order  differences  information  can  also 
be  appended. 

'I'be  methods  presented  have  also  been  investigated 
for  the  case  of  a  diffract  ion- limited  noncoherent  imag¬ 
ing  system.  The  discrete  formulation  of  the  nonco¬ 
herent  case  possesses  an  effectively  higher  rank  than  the 
coherent  case  as  a  result  of  the  larger  spatial  bandwidth. 
Thus  a  better  two-point  resolution  is  obtained  as  shown 
in  Fig.  (i  where  t  he  restored  image  (chain-  dot  curve)  is 
almost  identical  to  the  original  image  (solid  line),  and 
the  dotted  curve  represents  the  measured  DL  image. 
'I'be  discrete  spectrum,  similar  to  the  continuous 
transfer  function,  takes  on  intermediate  values  between 
zero  and  unity.  Therefore,  the  elimination  of  noisy 
image  components  in  the  presmoothing  step  is  not  as 
effective.  Since  the  image  components  retained  will 
represent  much  of  the  measurement  noise,  this  is 
demonstrated  in  Fig.  7  where  the  noncoherent  DL  image 
is  shown  by  the  dotted  curve.  'I'be  restored  image 


Kip.  (i.  Noncoherent  rest  oral  inn  of  two  impulses  separated  by  one- 
tenth  ot  the  Kayleiph  distance. 


Kip.  7.  Noncoherent  restoration  of  1  wo  impulses  separated  by  one- 
half  of  the  Ravleiph  distance  with  SNR  of‘JI  dB. 
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(chain-dot)  does  not  approximate  the  original  image 
(solid  line)  as  well  as  for  the  coherent  case.  Although 
the  SNR  is  I  dB  larger  than  in  that  of  the  coherent  case 
(Fig.  .r>),  the  restoration,  although  acceptable  for  many 
applications,  is  not  as  close  as  in  the  corresponding  co¬ 
herent  example. 


III.  Summary  and  Conclusions 

A  new  approach  to  image  restoration  has  been  in¬ 
troduced.  The  method  is  highly  immune  to  measure¬ 
ment  errors  and  noise.  The  ill-conditioned  nature  of 
the  problem  is  accounted  for  in  such  a  way  as  to  reduce 
the  amount  of  computations  necessary.  The  capability 
of  the  method  to  yield  superresolving  restorations  is 
demonstrated.  This  property  of  increasing  the  reso¬ 
lution  of  a  measured  image  beyond  that  dictated  by  the 
diffraction  limit  is  demonstrated  even  in  the  presence 
of  very  significant  noise.  The  method  is  particularly 
well-suited  for  impulsive  image  applications. 

This  paper  is  based  on  a  portion  of  a  dissertation 
submitted  by  R.  -I.  Mammone  in  partial  fulfillment  of 
the  requirements  for  the  I’h.I).  in  Electrical  Engineering 
to  the  City  University  of  New  York,  N.Y. 

This  work  was  supported  in  part  by  a  grant  under 
AFOSR  77-3217.  Some  of  this  material  was  presented 
at  the  1980  OSA  Annual  Meeting,  Chicago. 
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ABSTRACT 


A  method  of  restoring  the  discrete  Fourier  transform  (DFT) 
spectrum  of  a  diffraction-limited  (DL)  image  from  a  narrow 
observation  segment  of  the  DL  image  is  presented.  The  DL 
spectral  restoration  process  is  the  dual  of  the  more  common  DL 
image  restoration  process  with  the  roles  of  the  frequency  and 
space  reversed.  Applications  of  a  spectrum  restoration  include 
increasing  the  field  of  view  of  existing  imaging  systems  and 
extracting  precise  frequency  components  of  a  large  DL  image 
using  only  a  small  segment  of  the  entire  image.  This  method 
could  also  be  employed  for  image  data  compression  which  is  of 
interest  in  digital  video  applications.  Several  differences 
between  the  implementations  of  the  image  and  the  spectrum 
restoration  processes  are  described.  The  estimate  is  constrained 
to  have  an  upper  bound  on  the  number  of  frequency  components 
contained  in  the  Fourier  spectrum.  The  bound  is  the  number  of 
samples  acquired  at  the  Nyquist  rate  for  the  length  of  the 
image.  The  magnitude  of  the  discrete  Fourier  transform  (DFT)  < 

Bpeqtfum  la  also  bounded.  These  constraints  define  a  large 
number  of  possible  solutions.  The  desired  solution  is  then 


selected  such  that  the  distance,  defined  in  a  function-theoretic 
sense,  between  the  measured  and  the  estimated  image  is  an 
optimum.  A  number  of  such  measures  are  investigated.  Numerical 
experiments  show  that  this  approach  yields  results  that  are 
highly  immune  to  measurement  noise. 

*  This  paper  is  based  on  a  portion  of  a  dissertation  submitted 
by  R.J.Mammone  in  partial  fulfillment  of  the  requirements  for 
the  Ph.D.  degree  in  Engineering  to  the  City  University  of  New 
York.  This  work  was  supported  in  part  under  a  grant  from  the  Air 
Force  Office  of  Scientific  Research  AFOSR  81-0169  and  a  contract 
from  the  Rome  Air  Development  Center  F19628-80-C-0095 

**  presently  with  the  Department  of  Electrical  Engineering, 
Rutgers  University,  Piscataway,  N.J.  08854 
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INTRODUCTION 


The  finite  aperture  of  any  physical  imaging  system  eliminates 
the  high  spatial  frequency  components  of  the  object  from 
appearing  in  the  image.  The  lack  of  high  frequency  detail 
results  in  a  loss  of  resolution  in  the  observed  image.  It  lias 
been  shown  [1],  that  for  an  object  of  finite  extent,  an  exact 
restoration  of  the  object  from  the  d i f f r act  ion- 1 i mi  ted  (DL) 
image  is  possible.  The  restoration  of  the  DL  object  can  be 
viewed  as  a  continuation  of  the  spatial  Fourier  spectrum  beyond 
the  spatial  cutoff  frequency  imposed  by  the  DL  system.  However, 
numerical  implementation  of  DL  image  restoration  process  are 
highly  unstable  in  the  presence  of  measurement  noise.  The 
restoration  process  can  be  stabilized  by  imposing  additional 
constraints  [2].  Many  restoration  methods  are  available  such  as 
those  of  Frieden  [2],  Schell  [3],  Biraud  [4],  Jansson  [51  and 
Burg  [6].  Recently,  Howard  [7]  has  demonstrated  a 
computationally  inexpensive  method  using  a  least-squares 
approach.  This  approach  was  further  elaborated  on  by  Rushforth 
et  al.[8].  A  new  method  of  obtaining  an  increase  of  image 
resolution  using  linear  programmming  techniques  has  recently 
been  given  by  Mammone  and  Eichmann  [9]. 

In  this  paper,  the  dual  of  the  DL  image  restoration  problem  is 
considered.  Here  the  extrapolation  of  a  finite  segment  of  the  DL 


54) 


(i.e.  spatial  bandlimited)  image  data  in  the  presence  of 
measurement  noise  is  presented.  The  restoration  of  the  discrete 
Fourier  transform  (DPT)  spectrum  may  be  interpreted  as  an 
extrapolation  of  the  truncated  spatial  image.  Application  of  a 
spectrum  restoration  include  increasing  the  field  of  view  of 
existing  imaging  systems  and  extracting  exact  frequency 
components  of  a  large  DL  image  when  only  a  smaller  image  segment 
of  the  entire  image  is  available.  This  reduction  would  be  also 
of  interest  in  image  data  compression  applications  such  as  the 
transmission  of  digital  video  images.  While  spectral  restoration 
yields  the  same  mathematical  formulation  as  that  obtained  for  DL 
image  restoration,  with  the  role  of  space  and  spatial  frequency 
reversed,  there  are  significant  differences  as  well.  Here  the 
unknowns  are  Fourier  coefficients  which  are  complex  numbers. 

This  fact  appears  to  double  the  numerical  complexity  of  the 
problem.  Also,  the  circulant  convolutional  operator  which 
results  from  approximating  the  exact  Fourier  transform  with  a 
DFT  leads  to  additional  considerations.  Finally,  while  the 
restored  DL  object  is  non-negative,  the  restored  DFT  spectrum 
need  not  be  real  nor  non-negative.  Thus,  the  powerful 
non-negativity  constraint  used  in  the  DL  image  restoration  con 
not  be  directly  used. 

The  motivation  of  this  research  is  to  provide  high  resolution  t 

frequency  estimates  of  data  available  only  on  an  incomplete 
observation  lenqth  using  the  computationally  efficent  FFT 

i 


algorithms.  In  order  to  simplify  the  discussion  1-D  images  are 
considered.  For  a  DL  (band  1 imi ted )  image  truncated  to  lie  within 
an  interval  !J,  the  spatial  frequency  resolution  is  limited  by 
the  uncertainty  principle  [12].  The  lack  of  spectral  resolution 
is  related  to  the  infinite  limits  of  integration  in  the 
definition  of  the  Fourier  transform,  Since  the  image  is 
available  over  a  finite  observation  length,  the  calculated 
Fourier  transform  is  distorted.  It  is  convolved  with  a  sine 
function  with  mainlobe  width  D£ .  For  the  case  of  two  different 
sinusoidal  images  the  uncertainty  principle  states  that  the 
difference  in  the  frequency  between  the  images  can  not  be  less 
than  D£  otherwise  there  will  be  a  significant  overlap  of  the 
transforms.  This  overlap  would  not  permit  the  two  separate 
responses  to  be  distinguishable.  Thus,  the  frequency  resolution 
is  1 imi ted  to 

D£  ^  1/S  (1) 

In  many  practical  situations,  the  image  is  not  available  for  an 
interval  of  sufficient  length  S  for  the  desired  frequency 
resolution  D£.  It  is  then  desirable  to  process  the  spectrum  of 
the  observed  signal  such  that  the  spectral  resolution  is  greater 
than  that  given  by  Cq.(l).  If  there  is  no  source  of  measurement 
error,  it  is  well  known  [1]  that  the  spectrum  can  exactly  be 
restored.  However,  there  is  always  noise  and/or  error  present 
due  to,  if  nothing  else,  the  finite  numerical  precision 
necessary  to  carry  out  the  computation.  Therefore,  noise  or 
measurement  error  must  be  considered  as  part  of  the  formulation 
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of  the  problem.  In  this  case,  the  restoration  process  becomes 
highly  unstable. 

It  has  been  [2]  found  that  the  imposition  of  a  priori 
constraints,  such  as  nonnegativity  of  the  estimate,  will 
stabilize  the  restoration  process.  Since  the  maximum  number  of 
independent  equations  that  can  be  generated  is  equal  to  the 
number  of  spatial  samples  r  acquired  at  the  Nyquist  rate,  we 
shall  restrict  the  spectrum  restoration  method  co  situations 
where  the  number  of  frequency  components  is  less  then  or  equal 
to  r.  This  limit  on  the  dimension  of  the  estimate  has  been 
obtained  by  several  other  studies  [12,13,14],  The  sinusoidal 
frequency  components  will  be  further  assumed  to  be  harmonics  of 
the  same  fundamental  frequency.  This  latter  assumption  is  not 
necessary  but  is  made  to  facilitate  the  use  of  the  fast  Fourier 
transform  (FFT)  approximation  of  the  Fourier  spectrum.  The 
periodic  frequency  spectrum  is  obtained  by  zero-padding  the 
truncated  space  sequence  such  that  the  total  number  of  elements 
is  L,  an  integer  that  is  a  power  of  two.  This  step  allows  a 
finer  frequency  scale  to  be  generated. 

The  optimal  data  fitting  employs  linear  programming  (LP) 
techniques.  The  LP  method  provides  both  cost  and  time  effective 
ways  of  selecting  the  optimal  r-tuple  estimate.  LP  methods  have 
been  used  previously  for  the  general  image  restoration  problem 
[9,10,11].  The  advantage  of  the  LP  techniques  for  the  unstable 
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DL  image  restoration  problem  has  also  been  addressed  [9).  In 
this  paper  we  demonstrate  the  advantage  of  this  approach  for  the 
rank  deficient  spectrum  restoration  problem.  Since  the  rank,  the 
number  of  independent  equations,  is  less  then  the  number  of 
unknowns,  many  solutions  exist.  The  deficiency  of  the  rank 
occurs  due  to  the  DFT  low-pass  filter  effect  that  eliminates  all 
but  r  discrete  frequency  components.  This  efffect  is  in  contrast 
to  the  linear  discrete  convolution  case  [9]  which  yields  a 
system  of  equations  which  has  a  full  rank.  However,  because  the 
adjacent  equations  are  almost  identical,  when  a  finite 
arithmetic  machine  representation  is  used,  the  system  of 
equations  is  ill-posed,  i.e.  nearly  rank  deficient. 

FORMULATION 


For  the  following  discussion,  it  will  be  helpful  to  define 

negative  spatial  and  frequency  samples.  The  first  (L/2+1) 

elements  of  the  sequence  f  or  F,  correspond  to  the  positive 

n  k 

spatial  or  frquency  samples,  respectively,  in  ascending  order 
and  the  remaining  (L/2-1)  elements  correspond  to  the  negative 
spatial  or  frequency  samples,  respectively,  in  descending  order. 
This  is  consistent  with  the  idea  of  extending  the  sequences  in  a 
periodic  manner.  The  spatial  truncation  operation,  the  model  for 
a  short  observation  segment,  is  characterized  by  multiplication 
of  the  spatial  sequence  by  the  unit  rectangle  sequence  Rr  where 
r  Is  the  number  of  adjacent  unity  elements  centered  about  the 
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zero  element  and  with  all  other  elements  are  equal  to  zero.  The 

DFT  of  the  rectangle  sequence  is  real  and  even,  since  R^  is  real 

and  even,  and  is  similar  to  a  sampled  sine  function.  Let  the 

column  vector  h  be  defined  as  a  vector  whose  elements  are 

h  =  DFT {  R  }  (2) 

n  r 

We  shall  use  the  circulant  convolution  theorem  [15]  to  examine 
the  relationship  betv/een  the  truncated  and  the  exact  sequence. 
The  circulant  convolution  matrix  for  the  discrete  spatial 
truncation  operator  R  is  given  by  the  circulant  matrix 


and  therefore 

G  =  HF  +  N  (4) 

where  G  and  F  are  the  DFT  of  the  measured  and  the  actual 
sequence  vectors,  respectively,  and  N  is  the  complex  error 
(noise)  vector.  It  is  well  known  (15)  that  the  eigen-values  of  a 
circulant  matrix  M  are  the  coefficients  of  the  DFT  of  the  first 
row  of  the  circulant  matrix  and  the  eigen-vectors  are  the  DFT 
basis  vectors.  Therefore,  the  rank,  which  is  equivalent  to  the 
number  of  nonzero  eigenvalues,  of  I!  is  equal  to  the  number  of 
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non-zero  elements  r  in  the  rectangle  sequence  R  .  Since  the 
vectors  G ,  F  and  N  are  complex,  Eq.(4)  represents  2L  real  scalar 
equations.  But  all  the  spatial  sequences  are  real  and  therefore, 
the  DFT  sequences  must  have  symmetry.  The  symmetry  in  turn 
introduces  redundancies  in  Eq . ( 4 ) .  The  elimination  of  this 
redundancy  reduces  the  2L  to  L  real  equations.  In  order  to 
obtain  these  L  equations  we  must  examine  some  properties  of  the 
DFT. 

Similar  to  the  properties  of  the  Fourier  transform,  the  DFT  of  a 
real  sequence  possesses  an  even  real  and  an  odd  imaginary  part. 
Since  we  wish  to  use  the  computationally  efficient  radix  two  FFT 
algorithm,  the  sequences  must  consist  of  an  even  number  o L 
points.  The  existing  symmetry  can  be  displayed  in  the  following 
way.  There  is  one  element  for  the  dc,  the  zeroth,  component, 
(L/2-1)  negative  as  well  as  positive  elements  and  one  remaining, 
the  L/2  ,  element.  This  can  be  seen  by  noticing  that  the  L/2th 
row  of  the  DFT  matrix  is  a  sequence  of  alternating  positive  and 
negative  unities.  Thus,  the  center  and  the  dc  frequency 
components  are  always  real.  The  remaining  (L-2)  elements  consist 
of  complex  conjugate  pairs,  between  the  first  and  the  last 
elements  and  the  second  and  the  sr  • ->nd  to  the  last  elements  and 
so  on.  There  are  (L/2+1)  distinct  real  and  (L/2-1)  imaginary 
elements.  We  form  the  concatenated  sequence  of  the  distinct 
elements.  Let  the  real  and  the  imaginary  parts  of  the  ith 
element  be  denoted  by  the  subscript  Ri  and  Ii,  respectively, 
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then 


=  1  GR0 '  GR1 '  "  • 


GR ( L/2  +  1 ) '  Gr (L/2-1 ) '  * 


'  GI1  1 


=  1  FR0'  FR1'  **  '  FR ( L/2+1 ) '  FI (L/2-1 ) 


,  ..  ,  Fn  }T  (5) 


G 
F 

N  =  {  NR0'  "Rl» 
where  the  superscript  T  stands  Cor  the  transpose  operation.  The 


N, 


•  '  NR(r./2tl|'  NI(L/2-))' 


nh  1 


relation  Eq.(4)  can  now  be  written  in  reduced  form 

G  =  H  F  +  N  (6) 

where  if  we  decompose  the  degradation  matrix 


ll  = 


(7) 


I  e!G|P  I  °| 

where  A  is  a  ( L/2+1 ) x ( L/2+1 ) ,  B  is  a  ( L/2+ 1 ) x ( L/2 -1 ) ,  C  and  D  are 
( L/2-1 ) x ( L/2-1 )  submatrices  of  M  and  e  and  p  are  unused  column 
vectors.  If  we  denote  the  null  vector  and  matrix  0  and  0, 
respectively,  and  define 


B  =  {  0  B  0  } 


(8) 


C  =  {  0  0  C  }  (9) 

D  =  {  0  0  D  }  (10) 


then 


II 


A  +  B  j  0 
I 

_ l_. 

I 


(ID 


Thus  Eq.(6)  represents  L  real  equations  in  the  2L  variables  F  and 

A  A 

N  .  Here,  H  denotes  the  space  truncation  operator  in  the  DFT  , 
domain.  This  underdetermined  system  of  equations  has  many 
solutions.  For  example,  the  pseudo-inverse  solution  might  be  used 


i 
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(16),  although  this  would  produce  an  unconstrained  estimate.  In 
order  that  we  may  be  able  to  impose  constraints  on  and 
simultaneously  choose  an  optimal  solution,  that  is  one  close  in 
some  sense  to  the  measured  image,  we  use  a  method  of  optimization 
known  as  linear  programming. 


CONSTRAINED  OPTIMAL  SOLUTIONS 

The  basic  linear  programming  (LP)  problem  is  formulated  in  the 
following  way  [17]  : 

T 

Minimize  the  cost  function  :  c  x 

Subject  to  :  Ax  ^  h  (12) 

x  -  0 

where  c  and  x  are  the  known  cost  and  unknown  M-d imens iona 1 
variable  vectors,  respectively,  b  is  an  r-d imensional  constraint 
vector  and  A  is  a  rxM  constraint  matrix.  The  most  often  used  LP 
technique  is  the  simplex  method  [17].  The  output  obtained  from  a 
simplex  LP  algorithm  will  correspond  to  one  of  the  following  three 
situations:  either  there  is  no  feasible  solution,  that  is,  all  of 
the  constraints  cannot  simultaneously  be  satisfied,  or  the  optimal 
solution  is  feasible  but  is  unbounded.  The  third  alternative  is 
that  the  solution  is  feasible  and  bounded  and  then  an  optimal 
solution  is  provided. 

i 

The  first  step  in  the  simplex  method  is  to  convert  the  inequality 
constraints  into  equalities.  This  step  is  accomplished  by 


introducing  additional  variables,  called  slack  variables.  The 
resulting  system  of  equations  will  have  more  unknowns  than 
equation  leading  to  an  underdetermined  system  of  equations.  From 
the  many  possible  solutions  the  simplex  method  seeks  a  solution 
that  minimizes  the  cost  function.  The  basic  promise  of  LP  is  that 
the  solution  vecor  x  will  consist  of  at  most  r  non-zero  elements, 
where  r  equals  to  the  number  of  independent  equations,  and  the 
remaining  M-r  elements  will  be  equal  to  zero.  In  a  vector  space 
interpretation  of  LP,  the  inequalities  define  intersecting  hyper- 
planes  which  form  a  region  of  possible  solutions.  This  region  is 
an  r  dimensional  polygon  or  simplex.  The  simplex  will  always  lie 
in  the  positive  octant  of  the  vector  space.  The  cost  function  also 
defines  a  hyper-plane  for  a  fixed  cost  value.  For  this  value  to  be 
minimum,  the  cost  function  hyper-plane  must  lie  on  an  edge  of  the 
simplex.  The  coordinates  of  each  edge  contain  at  most  r  numbers. 
The  simplex  method  is  an  iterative  technique  which  begins  with  a 
feasible  solution,  at  a  particular  edge  of  the  simplex  and, 
progresses  to  an  adjacent  edge  in  such  a  way  as  to  yield  a 
reduction  in  the  cost  value  until  no  further  reduction  is 
possible.  The  coordinates  of  the  resulting  edge  yields  the  optimal 
solution. 

In  a  LP  formulation  of  the  spectral  estimation  problem,  we  now 
take  Eq.(fi)  to  to  be  part  of  the  constraint  matrix.  In  order  to 
guarantee  that  all  the  variables  are  non-negative,  a  necessary 
requirement  in  LP,  we  next  define  the  vectors  in  Eq. (6)  in  terms 


^  * 

of  the  nonnegative  component  vectors  F  ,  F  ,  N  and  N 


that 


F  =  F  -  F 


N  =  N  -  N 


If  we  examine  the  sum 

N  +  +  N  =|g  -  H  F  |  (15) 

we  observe  that  Eq  .(15)  is  the  1^  norm  of  the  error  between  the 
spectrum  of  the  measured  signal  G  and  the  spectrum  of  the 
spatially  truncated  image  estimate  F  .  This  is  the  error  measure 
to  be  minimized.  Thus  the  LP  formulation  of  the  1^  spectral 


estimation  problem  is 
Minimize  : 
Subject  to 


N  -F  N 


=  H  (F  +-F  )+N  +-N 


F  +  ^  F  + 


F  +,  F  ",  N  +,  N  Ss  0 


The  1^  norm  is  not  the  only  measure  of  error  that  can  be 

minimized.  The  1.  c  norm  or  maximal  deviation  can  also  be 

l  n  E 

employed.  This  can  be  seen  by  defining  a  scalar  E  to  be  the 
absolute  value  of  the  maximum  element  of  the  error  vector  given  by 
Eq.(lS).  Thus  the  LP  formulation  of  the  1.^  spectral  estimation 
problem  is 

Minimize  :  F. 


Subject  to 


G  ^  II  (F  -F  ~)+Ee 
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G  -  H  ( F  +-F  ”)-Ee 

F  F  +  (17) 

max 


F  +,  F  ”,  E  -  0 

where  e  is  an  L  dimensional  unit  vector, 
n 

The  1^  norm  or  the  sum  of  the  squares  of  the  error  can  also  be 
minimized  using  a  formulation  similar  to  that  of  the  1  ^  norm.  Here 
quadratic  programming  methods  are  used.  Quadratic  programming 
techniques  provide  an  estimate  with  at  most  r  positive  elements 
with  all  other  elements  as  zero.  In  general,  quadratic  programming 
techniques  use  computationally  efficient  LP  methods  but  on  larger 
augmented  matrices.  These  methods  are  well  known  and  many 
published  computer  routines  are  available  [18,19].  There  is 
usually  more  than  one  least-square  estimate  to  an  underdetermined 
system  of  equations  like  Eq.(4).  For  example,  the  pseudo-inverse 
solution  [15]  provides  an  unconstrained  leastsquare  estimate  v/ith 
lowest  1^  norm  as  well  as  the  minimum  square-error.  Howard  [7]  has 
demonstrated  a  constrained  least-squares  approach  for  the  dual 
problem  incorporating  a  penalty  function  to  force  the  solution  to 
become  nonnegative.  Rushforth  et  al.  [8]  developed  this  approach 
for  ill-posed  DL  image  restoration  problem. 

The  estimates  obtained  by  optimizing  the  1^,  l  and  l^n{.  norms  are* 
maximum  likelihood  estimates  when  the  noise  is  modeled  as  a  random 
variable  with  uniform,  gaussinn  or  exponential  probability  density 

f 

i 

l 
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functions  (pdf),  respectively  [20].  It:  was  found,  however,  that: 
for  the  low  noise  situation  of  interest,  the  pdf  of  t  lie  noise  did 
not  appreciably  effect  the  estimates.  In  terms  of  increased 
frequency  resolution,  the  1^  norm  consistently  offered  tin'  best 
estimate  i  r  respect,  ive  of  the  noise  pdf.  This  fact  can  be 
attributed  to  the  relatively  fewer  numerical  operations  needed  as 
well  as  the  strongest  bounds  on  the  individual  spectral  components 
for  the  1  ^  estimate.  The  constrained  minimal  ]  norm  using 
quadratic  programming  techniques  requires  much  more  numerical 
operations.  Also,  optimizing  the  1^  norm  is  equivalent  to 
minimizing  the  total  noise  spectral  energy.  While  this  is  a 
desirable  feature  in  general,  for  the  purpose  of  increased 
frequency  resolution,  it  tends  to  be  counterproductive.  The 
minimization  of  the  global  spectra)  error  noise  energy 
distribution  tends  to  overly  smooth  the  estimate  leading  to  a 
decrease  in  frequency  resolution. 


NUMERICAL  RESULTS 


The  observed  spectrum  is  a  smoothed  version  of  the  desired 
spectrum.  Fig.  1  illustrates  the  inadequacy  of  simply  taking  the 
FFT  of  a  truncated  sequence.  Mere  as  well  as  in  all  subsequent 
Figures  the  magnitude  of  the  Fourier  spectrum  is  plotted.  The 
dotted  line  represents  a  32  point  FFT  of  23  samples  of  a  cos  (w^x) 
sequence,  where  is  the  fundamental  frequency  2pi/32.  Mere,  only 
seven  data  samples  have  been  replaced  by  zeros.  The  solid  line 
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illustrates  the  exact  spectrum  of  the  full  32  sample  sequence.  The 
lack  of  resolution  conceals  the  fact  that  only  the  fundamental 
frequency  is  present.  The  size  of  the  FFT  as  well  as  the  number  of 
samples  acquired  are  known  a  priori.  This  information,  together 
with  the  bounds  on  the  amplitude  spectrum,  is  used  to  estimate  the 
actual  spectrum.  Fig.  2  illustrates  the  performance  of  the  1  ^  norm 
when  only  seven  out  of  the  thirty-two  samples  of  the  fundamental 
cosinusoid  is  available.  The  dotted  line  curve  represents  the 
overly  smooth  estimate  of  the  direct  FFT  spectrum.  The  minimal  1^ 
estimate,  the  dotted  line,  is  very  close  to  the  original  spectrum. 

The  dotted  line  curve  lies  completely  on  the  solid  line  curve. 

Thus,  a  complete  restoration  is  possible  even  when  only  seven  out 
of  thirty  two  samples  are  available.  Fig.  3  illustrates  the 
extrapolated  space  domain  sequence  co r respond i  vj  to  Fig.  2.  The 
acquired  samples  are  illustrated  by  the  solid  line  segment  and  the 
extrapolated  curve  by  the  dotted  line  curve. 

The  method  is  quite  robust  to  measurement  noise.  Fig.  4 
illustrates  the  restoration  of  the  fundamental  cosinusoidal 
waveform  from  fifteen  samples  of  a  thirty-two  point  FFT  with  white 
gaussian  noise  added.  The  variance  of  the  noise  is  one  tenth  of 
the  amplitude  of  the  function.  Thus,  an  increase  on  the  order  of  a 
factor  of  two  in  resolution  is  obtained  even  in  the  presence  of 
very  significant  measurement  noise.  The  corresponding  spatial  t 

extrapolation  is  depicted  in  Fig.  5.  Here  the  incomplete 
observation  is  shown  by  the  dotted  line  curve.  The  extrapolated 
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waveform  (the  chain  dotted  curve)  follows  the  exact  curve  (solid 
line)  very  closely.  In  many  applications  it  is  desired  to  resolve 
two  frequency  components  spaced  very  close  together.  Fig.  f> 
demonstrates  the  performance  of  the  1^  norm  in  this  case.  Mere  the 
first  and  the  fifth  harmonics  are  sampled  with  eleven  out  of  the 
possible  thirty-two  samples.  The  direct  application  of  the  FFT 
provides  the  blurred  spectrum  indicated  by  the  dotted  curve.  The 
estimated  Fourier  spectrum  is  represented  by  the  solid  line  where 
again,  the  estimate  follows  the  exact  spectrum  so  closely  that  the 
curves  are  indistinguishable. 

SUMMARY  AND  CONCLUSIONS 

A  new  method  of  extrapolating  a  DL  image  from  a  measured  spatially 

truncated  image  is  presented.  The  method  directly  addresses  the 

issue  of  the  finite  degree  oE  freedom  of  the  DL  image.  The 

feasible  image  estimates  are  constrained  to  have  an  upper  bound  on 

the  number  of  frequency  components  comprising  the  image.  A 

geometric  approach  has  been  taken,  that  is,  various  error  measures 

were  optimized  subject  to  given  constraints.  Three  error  norms 

were  discussed:  the  1.,  1 „  and  1.  c  norm.  This  method  allows  the 

1  2  inf 

addition  of  other  linear  constraints  on  the  desired  image.  That 
is,  any  inequality  possessing  a  linear  combination  of  the  unknown 
samples  on  one  side  of  the  inequality  and  a  scalar  value  on  the 
other  side  can  also  be  ar.  .  The  inclusion  of  constraints  of  this 
type  have  been  found  to  stabilize  the  otherwise  unstable  problem. 
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Numerical  results  demonstrate  the  increase  in  resolution  as  well 
as  its  robustness  to  measurement  noise.  In  terms  of  resolution, 
the  1^  norm  was  found  to  give  consistently  the  best  spectral 
estimates. 
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Figure  1  The  32  point  FFT  of  23  sample  points  of  cos(w.x)  (dotted 

Li 

line  curve)  and  the  exact  DFT  spectrum  (solid  line) 

Figure  2  The  restoration  of  the  DFT  of  cos  (w^x)  from  the  32  point 
FFT  of  7  spatial  domain  sample  points  (dotted  line  curve  )  and  the 
estimated  spectrum  (solid  line) 

Figure  3  Ex t r apo 1  a t i on  of  cos(w^x)  in  the  spatial  domain  (dotted 
line  curve)  from  7  sample  points  (solid  line) 

Figure  A  Restoration  of  the  cos(wrx)  from  15  spatial  sample  points 

Lj 

(dotted  line  curve  )  with  white  gaussian  noise.  The  noise  variance 
is  0.1.  The  restored  spectrum  is  depicted  by  a  solid  line. 

Figure  5  Extrapolation  of  cosfw^x)  (chain-dot  curve)  in  the 
spati  1  domain  corresponding  to  frequency  domain  plot  of  Fig. A. 

The  dotted  curve  indicates  the  incomplete  observation  with 
measurement  error.  The  extrapolated  curve  (chain-dot)  follows  the 
actual  waveform  closely  (solid  line  curve)  very  closely 
Figure  6  Restoration  of  cos  (wrx)  +  cos  (5w,x)  using  a  32  point 

U  L 

FFT  from  11  spatial  domain  sample  points.  The  DFT  is  depicted  with 
dotted  and  the  restored  waveform  is  shown  with  the  solid  curve. 


Figure 


78) 


292-09 


White-light  prefiltering  for  real-time  digital  image  transmission  of  still  and  moving  color,  video  iiroges 
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Abstract 

Sign  1 1  -dependent  noise  rncri  into  rod  when  samiding  video  signals  at  lev  sampling  rate,  will  le  using  an  adap¬ 
tive  video  del  La  modulator  Imth  .is  a  source  encoder  .and  a  video  digitizer,  is  cunhiLcd  by  using  a  white-light 
reflective  transform  opt  ical  preprocessor.  Whi  le  the  white-light  preprocessor  works  on  black  and  white  images, 
its  main  advantage  is  that  it  can  simultaneously  process,  using  a  single  source,  many  color  channels.  Further, 
the  relative  lack  of  both  spatial  and  temporal  coherence  aids  in  reducing  speckle  and  interference  noise 
effects.  Experimental  results  on  color  video  itrages,  that  have  been  preprocessed  using  a  white-light  optical 
transform  preprocessor  and  digitally  encoded  by  sampling  with  an  adaptive  delta  modulator  close  to  the  Nyquist 
rate,  are  presented. 


I nt  roduc  t  Lon 


Along  with  the  deve 1 0|xnent  of  sophisticated  monitoring  and  comuuni cation  equipment  the  need  for  real-time 
digital  video  transmission  is  apparent.  Digital  signal  transmission  offers  many  advantages  over  analog  .trans¬ 
mission  when  operating  over  a  channel  that  is  corrupted  by  noise.  This  is  due  to  the  fact  that  quantized 
levels  are  more  easily  di scriminated  from  the  effects  of  noise.  In  addition,  the  digital  format  allcvs  for 
time  multiplexing,  error-coding  and  correcting  as  well  as  encrypting,  procedures  that  are  not  available  on  an 
analog  channel.  However,  when  the  standard  four  MHz  video  signal  is  sampled  at  the  Nyquist  rate  and  encoded 
using  pulse  code  modulation  <TfM)  to  four  to  six  quantization  levels,  the  data  rate  becomes  48  -  64  Mbs.  In 
the  case  of  color  video  signals,  it  is  not  uncormv.n  to  require  85  -  90  Mias  digital  bandwidth  for  celor  PCM  en¬ 
coded  digital  video  signals.  This  high  digital  bandwidth  is  usually  unacceptable  in  many  applications. 

Various  digital  encoding  alternatives  have  been  explored  to  reduce  the  dicital  bandwidth  to  an  acceptable 
level.  One  such  solution  is  a  hard-wired  adaptive  delta  modulator  (ADM)  that  works  at  video  rates.  In  order 
to  compress  the  digital  hindwidth  to  its  ultimate  limit  it  is  necessary  to  operate  the  Aijf-1  close  to  its 
Nyquist  rate.  However,  when  an  ATM  runs  near  its  Nyquist  rate,  a  new  noise  source,  ami  effect  that  is  signal 
dependent  and  it  termed  edge  busyness,  appears.  Ihi s  noise  source  is  due  to  the  adaptive  property  of  the  A'M 
leading  to  oscillator/  sampling  behaviour  near  the  Nyquist  rate.  For  an  ADM  this  noise  source  is  the  fundamen¬ 
tal  limitation  in  reducing  the  digital  video  bandwidth.  This  signal  dependent  noise  source  is  combated  by 
preprocessing  the  color  video  imacos  by  a  white-light  reflective  transform  analog  preprocessor.  While  the 
uhite-1  i gilt  preprocessor  works  on  black  and  white  images,  its  main  advantage  is  that  it  can  simultaneously 
process,  using  a  single  source.  minv  color  channels.  Further,  the  relative  lack  of  both  spatial  and  temporal 
coherence  aids  in  reducing  optical  speckle  and  interference  noise  effects.  Experimental  results  on  color 
a'Heo  images,  that  have  been  preprocessed  using  a  white-light  optical  transform  preprocessor  and  digitally 
encoded  by  sampling  writh  an  ADM  closr  to  the  Nyquist  rate,  is  presented. 


background 


In  an  ADM,  the  step  size  algorithm  is  signal  dependent .  The  algoritlim  of  the  ADM  used  in  cur  work  was 
developed  by  Schilling  and  Song  [l).  In  this  algorithm,  the  step  size  is  dependent  on  the  previous  trans¬ 
mitted  bit  and  is  given  as 
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where  E.  ,  is  the  transmitted  digit,  S.  ,  is  the  present  sample  of  the  input  to  the  encoder,  Y,  .  is  the  step 
size  of  tlac  AIM,  X.  .  is  the  encoder's  estimate  of  the  input  sample  and  Y  .  and  Y  is  the  minimiu  and 
maximum  allowable  step  size  of  the  ADM,  respectively.  Fig.  1  shews  the  rrPiationshlf^lWrng  the  clock  pulse, 
the  inpu*  signal  ,  the  estimate  arnl  the  output  hit  stream.  Notice  that  the  system  responds  to  the  rising  edge 
of  each  clock  pulse.  Observe  that  when  (he  esllmite  X,  .  is  less  than  the  sample  of  the  Input  signal  S,  ,, 
the  transmitter)  bit  is  a  one  aul  the  step  size  is  increased  by  the  factor  of  one  .and  a  ha  1  f .  TH>.S<  thek£stl- 
mate  rises  exponentially,  When  an  overshoot  occurs,  the  transmitted  bit  is  a  minus  one  and  the  step  size  de¬ 
creases  by  a  factor  of  one-half.  This  form  of  adapt tve  encoding  takes  the  human  observer  vision  Into  account. 
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On.'  to  tin’  different  ii.iti.il  coixli  l  i  mis  ,  canni  v  from  line  to  line,  .is  well  ns  frame  to  frame,  a  vertical 
edge  eneoi  mtei  ed  in  the  imago  will  occur  at  tiiflei  -it  jnsitions  in  tin  scan.  As  the  samjd  ini;  rate  decreases 
the  error  between  the  steps,  and  e,peci.illy  near  a  vertical  <xlgo ,  increases.  Fie.  2  illustrate'  what  can 
occur  when  the  encoder  encounters  a  high  contrast  vertical  edge  for  two  difliient  initial  predictor  vilues. 

Mere  we  sec  the  dotted  trace  reaches  the  correct  value  in  five  slips,  as  op|*r.ed  to  a  previ  ns  or  a  following 
scan  line  trace,  shown  as  a  solid  line,  that  takes  seven  steps,  lliis  error  in  tile  AOM  is  called  the  slope 
overload  noise.  The  difference  is  t.ot  highly  visible  when  the  snarling  rate  Is  hi^i,  but  as  tlvu  rate  Is 
lowered,  its  visibility  increases.  Fig.  1  shows  the  combi  ait  inn  of  one -dimensional  slope  overload  error  and 
the  l wn-d i mens i on  1 1  line  to  line  difference  at  the  sane  high  contrast  hmndartes  which  causes  an  image  de¬ 
gradation  km  wn  as  "edge  busyness".  As  a  function  of  time,  at  low  sa/npl  ing  rates,  the  edge  seems  to  "wiggle" 
about  a  vertical  line.  This  noise  not  only  degrades  image  quality,  hut  also  leads  to  visual  fatigue.  Pie 
slope  overload  noise  can  Iv  reduced  by  sampling  at  a  faster  rate.  However,  if  we  wish  to  compress  the  band¬ 
width,  it  is  not  an  acceptable  solution. 

Uhite-1  ight  preprocessor 

Hie  purpose  of  the  white-light  preprccessor  is  to  manipulate  the  spatial  frequencies  of  the  image.  Since 
the  high  spatial  frequencies  are  due  to  the  high  contrast  edges  and  the  high  contrast  edges  leads  to  edge 
busyness ,  the  selective  filtering  of  the  high  spatial  frequencies  can  effect  the  signal  dependent  noise. 

While  an  all  digital  implementat ion  of  Imth  transform  processing  and  ADM  encoding  is  possible,  in  general,  di¬ 
gital  image  transforms  are  time  -construing.  Analog  transform  processors,  because  they  operate  on  the  full 
image  at  the  spoil  1  of  light,  are  natural  candidates  for  real-time  imago  processing  applications. 

Highly  coherent  analog  spatial  transform  processors  suffer  the  disadvantage  that  they  are  sensitive  to 
positional  tolerance  errors,  dirt  noise  on  the  lenses  as  well  as  problems  with  speckle  noise.  Furthermore, 
for  color  imago  processing,  a  mciN'r  of  color  channels  must  bo  created  to  simultaneously  process  a  color  image. 
By  decreasing  both  spatial  and  temporal  coherence  of  the  source,  some  of  these  objections  can  be  circumvented. 

It  is  known  that  a  white-light  source  can  actpii  re  spatial  coherence  by  passing  its  output  through  a  pinhole. 
In  this  case,  the  light  appears  to  emanate  from  a  single  point  source.  The  light  radiating  from  the  pinhole 
has  a  weaker  intensity,  as  compared  Co  a  laser  source,  anil  its  region  of  spatial  coherence  is  rather  restric¬ 
ted.  However,  using  a  low  light  level  camera  and  constricting  the  site  of  the  input  images,  these  problems 
can  be  overcome. 

The  arrangement  for  the  white-light  preprocessor  follows  the  physical  description  of  the  narrowband  trans¬ 
form  processor.  'Ihc  source  is  a  broad-hitd  2CO  W  Mercurv-Xonon  arc  lamp.  Its  output  is  first  spatially  fil¬ 
tered  by  focusing  it  on  a  pinhole.  The  size  of  the  pinhole,  370i.ni  diameter,  was  chosen  as  a  cixnprcxni  se .  A 
small  pinhole  would  block  mist  of  the  output  light  as  well  as  it  would  create  a  chromatic  diffraction  pattern. 

A  large  pinhole,  in  turn,  will  restrict  the  region  of  lateral  coherence.  Ihc  light  from  the  pinhole  is 
collimated  by  an  F/B  6"  diameter  lens.  The  angular  spread  of  the  collimated  beam  is  inversely  proportional  to 
the  lateral  coherence  of  the  white-light  system.  Instead  of  the  usual  lens  configuration,  we  have  chosen  a 
reflective  system.  There  are  two  possible  reflective  svstoms.  In  Fig.  Afai  a  fully  reflective  transform, 
while  in  Fig.  AO))  a  transmissive  filtering  system  is  shown.  While  the  second  configuration  leads  to  the  use 
of  off-axis  parabolic  mirrors,  we  have  used  the  first,  single  parabolic  mirror,  configuration.  For  this  geo¬ 
metry  the  object,  the  transform  and  the  image  planes  arc  all  formed  in  the  focal  plane  of  the  parabolic  mirror. 
By  mounting  the  object  off-axis  and  moving  it  out  of  the  focal  plane,  see  Fig.  5,  Che  three  regions  are  fully 
separated.  For  stationary  objects  the  color  transparency  was  innorsed  in  a  liquid-gate  holder.  For  moving 
images  a  16nm  pin-registered  film  projector  was  used.  The  projector  was  placed  on  its  vibration-isolated 
stand.  While  we  expected  problems  due  to  the  vibration  and  the  lack  of  index-matching  gate,  the  time-avera¬ 
ging  effect  of  the  moving  frame  tended  to  compensate  for  these  defects. 

The  spatial  filtering  is  perform'd  with  a  mirror  whoso  reflectivity  is  adjusted  by  a  transparency  in  front 
of  it.  Roth  continuous  as  well  as  binary  filtering  can  be  performed.  Binary  filtering  tends  to  degrade  image 
quality  abruptly.  Continuous-tone,  in  particular  lt'gari  thmic,  filtering  tends  to  work  better,  Hie  prepro¬ 
cessing  is  completed  by  converting  hack  the  conditioned  imige  no  1  focusing  it  on  a  color  video  camera.  The 
three  color  channels;  R,  0  aid  B  are  encoded  as  we] 1  as  decoded  by  three  hard-wired  video  delta  modulators. 

The  resultant,  images  can  then  he  viewed  either  on  a  color  monitor  or  recorded,  using  an  RGB  to  NTSC  col  -r 
signal  convertor,  on  a  video  recorder.  Stat  binary  and  noving  color  images,  using  8  Mbs  per  color  channel, 
can  N'  digitally  transmitted  with  good  fidelity. 


Kxperi  merit  a  1  rc'-nlts  were  presented  on  the  real-time  encoding  oi  moving  and  stationary  color  video  images. 
To  compress  the  digital  bandwidth,  a  white-light  analog  prepnxressor  has  been  utilised  to  condition  the  video 
Image  for  digital  cnccxling.  The  preprocessor  allows  the  enctxiing  of  color  vitlco  irmges  with  8  Mbs  per  color 
channel . 
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